text:
function y = lnCumGaussian(x)

% LNCUMGAUSSIAN log cumulative distribution for the normalised Gaussian.
% FORMAT
% DESC computes the logarithm of the cumulative Gaussian
% distribution.
% ARG X : input position.
% RETURN y : log probability of the value under the cumulative
% Gaussian.
%
% SEEALSO : erf, erfcx, cumGaussian, lnDiffCumGaussian, gaussOverDiffCumGaussian
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% NDLUTIL

index = find(x< 0);
if length(index)
  y(index) = -.5*x(index).*x(index) + log(.5) + log(erfcx(-sqrt(2)/2* ...
                                                    x(index)));
end
index = find(x>=0);
if length(index)
  y(index) = log(cumGaussian(x(index)));
end
y=reshape(y, size(x));

author: SheffieldML

repo: GPmat


text:
function [f] = psin(n,z)
%Psin   Arbitrary order Polygamma function valid in the entire complex plane.
%
%                         d^(n+1)
%        polygamma(n,z) = --------log(Gamma(z))
%                         dz^(n+1)
%
%usage: [f] = Psin(n,z)
%
%        if n is 0 or absent then f will be the Digamma function.
%        if n=1,2,3,4,5 etc then f will be
%        the tri-, tetra-, penta-, hexa-, hepta- etc gamma functon
%        Real(n) must be zero or positive.
%
%tested under versions 6.0 and 5.3.1
%
%        Z may be complex and of any size.
%
%        This program uses the partial fraction expansion of the
%        derivative of the Log of an excellent Lanczos series approximation
%        for the Gamma function. Accurate to about 12 digits.
%
%example: psin(101, -45.6-i*29.4)
%         is near 12.5 + 9*i
%
%example: psin(10, -11.5-i*0.577007813568142)
%         is near a root of the decagamma function
%
%example: x=[1:0.005:1.250]'; [x gamma(x) log(gamma(x)) psin(0,x) psin(1,x)]
%         recreates Table 6.1 page 267 from A&S
%
%example: x=[1:0.01:2.00]'; [x psin(2,x) psin(3,x)]
%         recreates Table 6.2 page 271 from A&S
%
%example: x=1; y=[0:0.1:10]'; f=psin(0,x+i*y); [y real(f) imag(f)]
%         recreates Table 6.8 page 288 from A&S
%
%example: x=2; y=[0:0.1:10]'; f=psin(0,x+i*y); [y real(f) imag(f)]
%         recreates the last part of Table 6.8 page 293 from A&S
%
%References: C. Lanczos, SIAM JNA  1, 1964. pp. 86-96
%            Y. Luke, "The Special ... approximations", 1969 pp. 29-31
%            Y. Luke, "Algorithms ... functions", 1977
%            J. Spouge,  SIAM JNA 31, 1994. pp. 931
%            W. Press,  "Numerical Recipes"
%            S. Chang, "Computation of special functions", 1996
%
%
%see also:   GAMMA GAMMALN GAMMAINC
%see also:   mhelp psi
%see also:   mhelp GAMMA

%Paul Godfrey
%pgodfrey@conexant.com
%July 22, 2004
%see gamma for calculation details...

%this routine still works even if n is complex with real(n)>=0
%don't know what the resulting function is called though

% we have bit of a problem for real(z)<0
% we could use the reflection formula, eq #6.4.7 in A&S
% but arbitrary order derivs of cot are cumberson to compute
% (it generates a polynomial in powers of cot)
% so instead we will use a modification of eq #6.4.6
% and call this program recursively to shift z by 500 until real(z)>=0

if nargin==1
   z=n;
   n=0;
end

if n==0
   f=psi(z);
   return
end

if real(n)<0
   error('Invalid Polygamma order')
end

sizeofz=size(z);

isneg=find(real(z)< 0);
isok =find(real(z)>=0);

negmethod=1;
if ~isempty(isneg)
   if negmethod==0
      zneg=z(isneg);
      gneg=psin(n,zneg+1); % recurse if to far to the left...
      hneg=-(-1).^n.*gamma(n+1).*zneg.^-(n+1);
      fneg=gneg+hneg;
   else
      zneg=z(isneg);
      %shift by, say, 500, to speed things up.
      m=500;
      gneg=psin(n,zneg+m);
      hneg=0;
      for k=m-1:-1:0
          hneg=hneg+(zneg+k).^-(n+1);
      end
      hneg=-(-1).^n.*gamma(n+1).*hneg;
      fneg=gneg+hneg;
   end
end
if ~isempty(isok)
   z=z(isok);
end

% the zeros of the Lanczos PFE series when g=607/128 are:

r=[ -4.1614709798720630-.14578107125196249*i;
    -4.1614709798720630+.14578107125196249*i;
    -4.3851935502539474-.19149326909941256*i;
    -4.3851935502539474+.19149326909941256*i;
    -4.0914355423005926;
    -5.0205261882982271;
    -5.9957952053472399;
    -7.0024851819328395;
    -7.9981186370233868;
    -9.0013449037361806;
    -9.9992157162305535;
   -11.0003314815563886;
   -11.9999115102434217;
   -13.0000110489923175587];

% the poles of the Lanczos PFE series are:
%p=[ 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -10 -11 -12 -13];

e=exp(1);
g=607/128; % best results when 4<=g<=5
h=1/2;

%compute tricky PFE expansion of the deriv of the log
%of the Lanczos series. The series Lanczos' coeffs manifest
%themselves in the zero locations. The residues are all +/- 1
%compare this to A&S page 259, eq# 6.3.16 and page 260, eq# 6.4.10

s=0;
for k=length(r)-1:-1:0
    s=s+(1./((z-r(k+1)).^(n+1))-1./((z+k).^(n+1)));
end
%what happens if n is not a positive integer?
s=(-1).^n.*gamma(n+1).*s;

zgh=z+(g-h);
if n==0
%  s=log(zgh)+(-g./zgh + s);
%  use existing more accurate digamma function if n=0
%  should never reach this code since we trapped it above
   f=psi(z);
   return
else
% do derivs of front end stuff
   s=(-1)^(n+1)*(gamma(n).*zgh.^(-n) + g*gamma(n+1).*zgh.^-(n+1))+s;
end

if ~isempty(isneg)
   f(isneg)=fneg;
end
if ~isempty(isok)
   f(isok)=s;
end

f=reshape(f,sizeofz);

return

%a demo of this program is

warning off
x=[-5:1/64:5]';

figure(1)
axis([min(x) max(x) -20 20])
grid on
hold on

y=[];
for n=0:6
    y(:,n+1)=psin(n,x);
end

plot(repmat(x,1,size(y,2)),y)

figure(2)
x=-10:1/16:10;
y=-5:1/16:5;
[X,Y]=meshgrid(x,y);
z=complex(X,Y);
f=psin(4,z);
g=log10(abs(f));
mesh(x,y,g)
rotate3d on

z=-76+54*i;
[z psin(98, z)]

disp('A zero of Psin1')
z = -0.412134547951937-i*0.597811942320597;
[z psin(1,z)]

warning on
return

% Include this complex psi function
% in case user doesn't have one

function [f] = psi(z)
%Psi     Psi (or Digamma) function valid in the entire complex plane.
%
%                 d
%        Psi(z) = --log(Gamma(z))
%                 dz
%
%usage: [f] = psi(z)
%
%tested under versions 6.0 and 5.3.1
%
%        Z may be complex and of any size.
%
%        This program uses the analytical derivative of the
%        Log of an excellent Lanczos series approximation
%        for the Gamma function.
%
%References: C. Lanczos, SIAM JNA  1, 1964. pp. 86-96
%            Y. Luke, "The Special ... approximations", 1969 pp. 29-31
%            Y. Luke, "Algorithms ... functions", 1977
%            J. Spouge,  SIAM JNA 31, 1994. pp. 931
%            W. Press,  "Numerical Recipes"
%            S. Chang, "Computation of special functions", 1996
%
%
%see also:   GAMMA GAMMALN GAMMAINC
%see also:   mhelp psi
%see also:   mhelp GAMMA

%Paul Godfrey
%pgodfrey@intersil.com
%July 13, 2001
%see gamma for calculation details...

siz = size(z);
z=z(:);
zz=z;

f = 0.*z; % reserve space in advance

%reflection point
p=find(real(z)<0.5);
if ~isempty(p)
   z(p)=1-z(p);
end

%Lanczos approximation for the complex plane

g=607/128; % best results when 4<=g<=5

c = [  0.99999999999999709182;
      57.156235665862923517;
     -59.597960355475491248;
      14.136097974741747174;
      -0.49191381609762019978;
        .33994649984811888699e-4;
        .46523628927048575665e-4;
       -.98374475304879564677e-4;
        .15808870322491248884e-3;
       -.21026444172410488319e-3;
        .21743961811521264320e-3;
       -.16431810653676389022e-3;
        .84418223983852743293e-4;
       -.26190838401581408670e-4;
        .36899182659531622704e-5];

n=0;
d=0;
for k=size(c,1):-1:2
    dz=1./(z+k-2);
    dd=c(k).*dz;
    d=d+dd;
    n=n-dd.*dz;
end
d=d+c(1);
gg=z+g-0.5;
%log is accurate to about 13 digits...

f = log(gg) + (n./d - g./gg) ;

if ~isempty(p)
   f(p) = f(p)-pi*cot(pi*zz(p));
end

p=find(round(zz)==zz & real(zz)<=0 & imag(zz)==0);
if ~isempty(p)
   f(p) = Inf;
end

f=reshape(f,siz);

return

%A demo of this routine is:
clc
clear all
close all
x=-4:1/16:4.5;
y=-4:1/16:4;
[X,Y]=meshgrid(x,y);
z=X+i*Y;
f=psi(z);
p=find(abs(f)>10);
f(p)=10;

mesh(x,y,abs(f),phase(f));
view([45 10]);
rotate3d;

figure(2);
ezplot psi;
grid on;

One=psi(2)-psi(1)
EulerGamma=-psi(1)

return

author: Sable

repo: mcbench-benchmarks


text:
classdef PTKAirwayRadiusApproximation < PTKPlugin
    % PTKAirwayRadiusApproximation. Plugin to approximate radius of branches in
    %     the airway tree
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     PTKAirwayRadiusApproximation uses a distance transform to estimate the
    %     radius of each branch in the airway tree. The results are shown as an
    %     output image.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Airway radius <BR>approximation'
        ToolTip = 'Shows airways coloured by their radius approximation'
        Category = 'Airways'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = true
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
    end

    methods (Static)
        function results = RunPlugin(dataset, reporting)
            results = dataset.GetTemplateImage(PTKContext.LungROI);
            airway_results = dataset.GetResult('PTKAirways');
            results = PTKAirwayRadiusApproximation.GetRadiusApproximationFromAirwayTree(airway_results.AirwayTree, results, reporting);
        end
    end

    methods (Static, Access = private)

        function results = GetRadiusApproximationFromAirwayTree(airway_tree, template, reporting)
            airway_segmented_image = PTKGetImageFromAirwayResults(airway_tree, template, false, reporting);
            dt_image = airway_segmented_image.RawImage;
            dt_image = dt_image == 0;
            dt_image = bwdist(dt_image);

            segments_to_do = airway_tree;
            min_voxel_size_mm = min(airway_segmented_image.VoxelSize);

            segmented_image = zeros(template.ImageSize, 'single');

            while ~isempty(segments_to_do)
                segment = segments_to_do(end);
                segments_to_do(end) = [];
                airway_points_local = template.GlobalToLocalIndices(segment.GetAllAirwayPoints);
                max_radius = max(dt_image(airway_points_local))*min_voxel_size_mm;
                segmented_image(airway_points_local) = max_radius;
                segments_to_do = [segments_to_do segment.Children];
            end

            results = airway_segmented_image.BlankCopy;
            results.ChangeRawImage(segmented_image);
            results.ImageType = PTKImageType.Scaled;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function s = frac_to_s ( f )

%*****************************************************************************80
%
%% FRAC_TO_S writes a positive fraction into a left justified character string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real F, the number to be written into the string.
%    F should be between 0.0 and 1.0.
%
%    Output, string S, a representation of F.
%
  if ( f < 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FRAC_TO_S - Fatal error!\n' );
    fprintf ( 1, '  The input fraction F was negative:\n', f );
    error (  'FRAC_TO_S - Fatal error!' );
  elseif ( 1.0 <= f )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FRAC_TO_S - Fatal error!\n' );
    fprintf ( 1, '  The input fraction F = %g was 1 or more:\n', f );
    error (  'FRAC_TO_S - Fatal error!' );
  end

  s2 = sprintf ( '%f', f );

  s = s2[2:];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef PTKImageType
    % PTKImageType. An enumeration used to specify how an image should be displayed
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    enumeration
        Grayscale,  % Image is a 16-bit integer greyscale image
        Colormap,   % Image is an 8-bit integer colormap, using the Lines colormap
        Scaled,     % Image is a single-colour floating point and will be scaled between its minimum and maximum values
        Quiver      % Image is a quiver plot; each voxel has 3 component values defining a 3D vector
    end

end


text:
% Fig. 6.21   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clear all;
%close all;
clf

num=1;
den=[1 2 1];
rlocus(num,den);
axis([-3 1 -2 2])
grid;
title('Fig. 6.21: Root locus of G(s)=1/(s+2)^2 vs K');


text:
% ASTROTIK by Francesco Santilli

function [orbits0,t0,mi,R,mi0,J] = checkscenario(scenario)

    if ~isa(scenario,'struct')
        error('Wrong class of an input argument.')
    end

    f(1) = isfield(scenario,'orbits0');
    f(2) = isfield(scenario,'t0');
    f(3) = isfield(scenario,'mi');
    f(4) = isfield(scenario,'R');
    f(5) = isfield(scenario,'mi0');

    if ~all(f)
        error('Missing field in scenario structure.')
    end

    orbits0 = scenario.orbits0;
    t0      = scenario.t0;
    mi      = scenario.mi;
    R       = scenario.R;
    mi0     = scenario.mi0;

    [J,D] = check(orbits0,2);
    check(t0,1)
    check(mi,1)
    check(R,1)
    check(mi0,0)

    if ~(D==3 || D==5)
        error('Wrong size of input arguments.')
    end

    if ~(length(t0)==J && length(mi)==J && length(R)==J)
        error('Wrong size of input arguments.')
    end

    if any(mi < 0)
        error('mi must be a positive value.')
    end

    if any(R < 0)
        error('R must be a positive value.')
    end

    if mi0 <= 0
        error('mi0 must be a strictly positive value.')
    end

end


text:
function varargout = QuadAnim4(varargin)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %{
 This GUI provides an "animation-like" view of simulation output.

    Additional quadcopter modeling and simulation materials available at:
    github.com/dch33/Quad-Sim

Copyright (C) 2014 D. Hartman, K. Landis, S. Moreno, J. Kim, M. Mehrer

License:
    This program 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.

    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 LESSER GENERAL PUBLIC LICENSE
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%{
Quadcopter "Animation" GUI.
(c) David Hartman, 2014

If you're looking in this code for "answers" or a quick fix...good luck.
Though I've tried to add some comments throughout, altering this code is
likely to be an unpleasant exercise.
If it's any consolation, it wasn't very fun to make this mess either :P
I hope to release a more clearly written and faster version of this
code sometime in the future. Or feel free to do it for me and email it to
me!! Please!?!
%}

%.

% QuadAnim4 MATLAB code for QuadAnim4.fig
%      QuadAnim4, by itself, creates a new QuadAnim4 or raises the existing
%      singleton*.
%
%      H = QuadAnim4 returns the handle to a new QuadAnim4 or the handle to
%      the existing singleton*.
%
%      QuadAnim4('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in QuadAnim4.M with the given input arguments.
%
%      QuadAnim4('Property','Value',...) creates a new QuadAnim4 or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before QuadAnim4_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to QuadAnim4_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help QuadAnim4

% Last Modified by GUIDE v2.5 30-Apr-2014 10:38:36

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @QuadAnim4_OpeningFcn, ...
                   'gui_OutputFcn',  @QuadAnim4_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before QuadAnim4 is made visible.
function QuadAnim4_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to QuadAnim4 (see VARARGIN)

% Choose default command line output for QuadAnim4
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);
handles.AZval = get(handles.AZslider,'Value')*360; % Get azimuth view in degrees
handles.ELval = get(handles.ELslider,'Value')*90;
handles.j     = 1; % Initialize animation frame counting variable
handles.skipFlag = 0;
% Move GUI to center of screen
movegui('center')
% This sets up the initial plot - only do when we are invisible
% so window can get raised using QuadAnim4.
A = evalin('base', 'yout');
% Generate the geometry used to draw the quadcopter
r = .5; d = 1.25; h = .25; %inches: rotor dia., quad motor distance from
% center of mass, and rotor height above arms (entirely cosmetic)
a = 1; b = 1; c = 0.2; % Top , bottom, and sides lengths of hub (?)

% Construct rotor representations
N = [d  0 h].';% m1 rotor center [X Y Z]
E = [0 -d h].';% m4 rotor center
W = [0  d h].';% m2 rotor center
S = [-d 0 h].';% m3 rotor center
Nr = circlePoints(N, r, 10); Nr = [Nr Nr(:,1)]; % Rotor blade circles
Er = circlePoints(E, r, 10); Er = [Er Er(:,1)];
Wr = circlePoints(W, r, 10); Wr = [Wr Wr(:,1)];
Sr = circlePoints(S, r, 10); Sr = [Sr Sr(:,1)];
% Motors connecting to center of blade circles
mN = [d,d;
      0,0;
      h,0];
mE = [0,0;
     -d,-d;
      h,0];
mW = [0,0;
      d,d;
      h,0];
mS = [-d,-d;
       0,0;
       h,0];
% Construct body plot points
bNS = [ d, -d;
        0,  0;
        0,  0]; %For drawing the body "X" shape
bEW = [ 0,  0;
        d, -d;
        0,  0];
% Body (HUB) Squares
Top = [ a/2,   0,-a/2,   0;
          0, b/2,   0,-b/2;
        c/2, c/2, c/2, c/2];
Bot = vertcat(Top(1:2,:),-Top(3,:)); % Bot is same as top just below the center of mass
NEB = [ a/2, a/2,   0,   0;
          0,   0, b/2, b/2;
        c/2,-c/2,-c/2, c/2]; % By the way... north east is actually north west from above since x is north and y is east :P
NWB = [ a/2, a/2,   0,   0;
          0,   0,-b/2,-b/2;
        c/2,-c/2,-c/2, c/2];
SEB = -NWB;
SWB = -NEB;

phi = A(1,4);
the = A(1,5);
psi = A(1,6);

% ROTATION MATRIX --- ZYX ROTATION (R = Rib)
R = [cos(psi)*cos(the) cos(psi)*sin(the)*sin(phi)-sin(psi)*cos(phi) cos(psi)*sin(the)*cos(phi)+sin(psi)*sin(phi);
       sin(psi)*cos(the) sin(psi)*sin(the)*sin(phi)+cos(psi)*cos(phi) sin(psi)*sin(the)*cos(phi)-cos(psi)*sin(phi);
       -sin(the)         cos(the)*sin(phi)                            cos(the)*cos(phi)];

% Rotate body frame velocity vector
U = A(:,7);
V = A(:,8);
W = A(:,9);
Vi = zeros(length(A),3);
MvMax= max(sqrt(U.^2+V.^2+W.^2));
Vb = 3/MvMax*[U(1), V(1), W(1)]'; % Scale velocity
Vi(1,:) = R*Vb; % Rotate velocity vector to inertial frame for plotting

% Support for X-configuration nifty trick (Requires that quadModel
% structure is in base workspace)
if (~evalin('base','quadModel.plusConfig'))
    Rz = [ sqrt(2)/2, sqrt(2)/2, 0;
          -sqrt(2)/2,sqrt(2)/2, 0;
                   0,          0, 1];
    Nr = Rz*Nr;
    Er = Rz*Er;
    Wr = Rz*Wr;
    Sr = Rz*Sr;
    mN = Rz*mN;
    mE = Rz*mE;
    mW = Rz*mW;
    mS = Rz*mS;
    bNS = Rz*bNS;
    bEW = Rz*bEW;
    Top = Rz*Top;
    Bot = Rz*Bot;
    NEB = Rz*NEB;
    NWB = Rz*NWB;
    SWB = Rz*SWB;
    SEB = Rz*SEB;
end

% Rotate body parts Via Initialized R
NrR = R*Nr;
ErR = R*Er;
WrR = R*Wr;
SrR = R*Sr;
mNr = R*mN;
mEr = R*mE;
mWr = R*mW;
mSr = R*mS;
bNSR = R*bNS;
bEWR = R*bEW;
TopR = R*Top;
BotR = R*Bot;
NEBR = R*NEB;
NWBR = R*NWB;
SWBR = R*SWB;
SEBR = R*SEB;

% Plot the box rotation and ang. velocity and inertial frame velocity
% vector
axes(handles.axes1)
plot3(bNSR(1,:),bNSR(2,:),bNSR(3,:),'b','LineWidth',3) % Body Arm
hold on
plot3(bEWR(1,:),bEWR(2,:),bEWR(3,:),'b','LineWidth',3) % Body Arm
plot3(mNr(1,:),mNr(2,:),mNr(3,:),'k','LineWidth',4) % Motor
plot3(mEr(1,:),mEr(2,:),mEr(3,:),'k','LineWidth',4) % Motor
plot3(mWr(1,:),mWr(2,:),mWr(3,:),'k','LineWidth',4) % Motor
plot3(mSr(1,:),mSr(2,:),mSr(3,:),'k','LineWidth',4) % Motor
plot3(NrR(1,:),NrR(2,:),NrR(3,:),'g') % Motor 1 blades
plot3(ErR(1,:),ErR(2,:),ErR(3,:),'k') % Motor 4 blades
plot3(WrR(1,:),WrR(2,:),WrR(3,:),'k') % Motor 2 blades
plot3(SrR(1,:),SrR(2,:),SrR(3,:),'g') % Motor 3 blades
grey = [0.5 0.5 0.5];
top = fill3(TopR(1,:),TopR(2,:),TopR(3,:),'r'); alpha(top,0.8); % Top Surface
bot = fill3(BotR(1,:),BotR(2,:),BotR(3,:),'g'); alpha(bot,0.8); % Bottom Surface
ne  = fill3(NEBR(1,:),NEBR(2,:),NEBR(3,:),'c'); alpha(ne,0.8); % North East surface
nw  = fill3(NWBR(1,:),NWBR(2,:),NWBR(3,:),grey); alpha(nw,0.8); % North West surface
sw  = fill3(SWBR(1,:),SWBR(2,:),SWBR(3,:),grey); alpha(sw,0.8); % South West surface
se  = fill3(SEBR(1,:),SEBR(2,:),SEBR(3,:),grey); alpha(se,0.8); % South East surface
axis square
xlabel('X')
ylabel('Y')
zlabel('Z')
xlim([-3 3])
ylim([-3 3])
zlim([-3 3])
view(handles.AZval,handles.ELval)
grid on

omi = zeros(length(A),3); % Initialize omega inertiaF points array (L(A)x3)
P = A(:,1); Q = A(:,2); Rw = A(:,3);
MombMax = max(sqrt(P.^2+Q.^2+Rw.^2)); % Calculate max magnitude of omb
omb = 3/MombMax*[P,Q,Rw].'; % Store and scale current omega bodyF
omi(1,:) = R*omb(:,1); % Rotate omegab to inertiaF store in omegai array
qp1 = quiver3(0,0,0,omi(1,1),omi(1,2),omi(1,3),'ro');
qp2 = quiver3(0,0,0,Vi(1,1),Vi(1,2),Vi(1,3),'k');
hold off
% legend([qp1 qp2],'Angular Velocity','I-Frame Velocity','Location',...
%      'SouthOutside')
%  pp1 = plot3(omi(1,1),omi(1,2),omi(1,3),'b.');
%  pp2 = plot3(Vi(1,1),Vi(1,2),Vi(1,3),'k.');
%  legend([qp1 qp2 pp1 pp2],'Angular Velocity','I-Frame Velocity',...
%      'Ang. Vel. Track','I-Frame Vel. Track','Location','SouthOutside')
% hold off

minX = min(A(:,10))*3.28;
minY = min(A(:,11))*3.28;
% minZ = min(A(:,12))*3.28;
maxX = max(A(:,10))*3.28;
maxY = max(A(:,11))*3.28;
maxZ = max(A(:,12))*3.28;

% Plot the three dimensional trajectory of the box
axes(handles.axes2)
X = A(1:1,10)*3.28; Y = A(1:1,11)*3.28; Z = A(1:1,12)*3.28;
scatter3(X,Y,Z,36,'blue')
hold on
fill3([minX-1 maxX+1 maxX+1 minX-1],...
    [minY-1 minY-1 maxY+1 maxY+1],...
    [0 0 0 0],'g');
%alpha(0.7);
xlabel('X (ft)')
ylabel('Y (ft)')
zlabel('Z (ft)')
xlim([minX-1 maxX+1])
ylim([minY-1 maxY+1])
zlim([-0.1 maxZ+1])
view(handles.AZval,handles.ELval)
axis square
grid on
hold off

% phi_A = A(:,4);
% the_A = A(:,5);
% psi_A = A(:,6);

% Choose default command line output for animationGUI2
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes QuadAnim4 wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = QuadAnim4_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in startStop.
function startStop_Callback(hObject, eventdata, handles)
% hObject    handle to startStop (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
plotAnim(hObject, handles) % Pass handle to start/stop button to plot function

% --- Executes on slider movement.
function ELslider_Callback(hObject, eventdata, handles)
% hObject    handle to ELslider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles.ELval = get(hObject,'Value')*90; % Set elevation view in degrees
handles.AZval = get(handles.AZslider,'Value')*360; % Set aximuth view in degrees
view(handles.axes1,handles.AZval,handles.ELval);
view(handles.axes2,handles.AZval,handles.ELval);
drawnow
guidata(hObject,handles)

% --- Executes during object creation, after setting all properties.
function ELslider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ELslider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function AZslider_Callback(hObject, eventdata, handles)
% hObject    handle to AZslider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles.AZval = get(hObject,'Value')*360; % Set azimuth view in degrees
handles.ELval = get(handles.ELslider,'Value')*90; % Set elevation view in degrees
view(handles.axes1,handles.AZval,handles.ELval);
view(handles.axes2,handles.AZval,handles.ELval);
drawnow
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function AZslider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to AZslider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes during object creation, after setting all properties.
function axes1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to axes1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate axes1

function points = circlePoints(center, radius, numberOfPoints)
% Helper function for plotting points
% Inspired by "circle()" from Peter Corke's MATLAB Robotics Toolbox
c = center.'; % [x;y] location of center
r = radius;
n = numberOfPoints;
% compute points around the circumference
th = (0:n-1)'/n*2*pi; % angles coresponding to each point
x = r*cos(th) + c(1); % x part
y = r*sin(th) + c(2); % y part
points = [x,y].';
    if length(c) > 2
        z = ones(size(x))*c(3); % z part
        points = [x, y, z].';
    end

function frameSkips_Callback(hObject, eventdata, handles)
% hObject    handle to frameSkips (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of frameSkips as text
%        str2double(get(hObject,'String')) returns contents of frameSkips as a double
handles.frameSkipVal = str2double(get(hObject,'String'));
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function frameSkips_CreateFcn(hObject, eventdata, handles)
% hObject    handle to frameSkips (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
handles.frameSkipVal = str2double(get(hObject,'String'));
guidata(hObject,handles);

function plotAnim(button, handles)
hObject = button;
mode = get(hObject,'String');
handles = guidata(gcbo);
j = handles.j;
if strcmp(mode,'Start')
    set(hObject,'String','Stop')
    guidata(hObject,handles)
    A = evalin('base', 'yout');
    tout = evalin('base', 'tout');
    frameSkipVal = str2double(get(handles.frameSkips,'String'))+1; % Size of steps to take for plotting animation (1 plots
    % every frame, 2 about every other, etc.

    r = .5; d = 1.25; h = .25; %inches: rotor dia., quad motor distance from
    % cm, and rotor height above arms (entirely cosmetic)
    a = 1; b = 1; c = 0.2;

    % Construct rotor representations
    N = [d  0 h].';% m1 rotor center
    E = [0 -d h].';% m4 rotor center
    W = [0  d h].';% m2 rotor center
    S = [-d 0 h].';% m3 rotor center
    Nr = circlePoints(N, r, 20); Nr = [Nr Nr(:,1)]; % Rotor blade circles
    Er = circlePoints(E, r, 20); Er = [Er Er(:,1)];
    Wr = circlePoints(W, r, 20); Wr = [Wr Wr(:,1)];
    Sr = circlePoints(S, r, 20); Sr = [Sr Sr(:,1)];
    % Motors connecting to center of blade circles
    mN = [d,d;
          0,0;
          h,0];
    mE = [0,0;
         -d,-d;
          h,0];
    mW = [0,0;
          d,d;
          h,0];
    mS = [-d,-d;
           0,0;
           h,0];
    % Construct body plot points
    bNS = [ d, -d;
            0,  0;
            0,  0]; %For drawing the body "X" shape
    bEW = [ 0,  0;
            d, -d;
            0,  0];
    % Body
    Top = [ a/2,   0,-a/2,   0;
              0, b/2,   0,-b/2;
            c/2, c/2, c/2, c/2];
    Bot = vertcat(Top(1:2,:),-Top(3,:));
    NEB = [ a/2, a/2,   0,   0;
          0,   0, b/2, b/2;
        c/2,-c/2,-c/2, c/2];
    NWB = [ a/2, a/2,   0,   0;
              0,   0,-b/2,-b/2;
            c/2,-c/2,-c/2, c/2];
    SEB = -NWB;
    SWB = -NEB;

    phi = A(1,4);
    the = A(1,5);
    psi = A(1,6);
    % Rib
    R = [cos(psi)*cos(the) cos(psi)*sin(the)*sin(phi)-sin(psi)*cos(phi) cos(psi)*sin(the)*cos(phi)+sin(psi)*sin(phi);
           sin(psi)*cos(the) sin(psi)*sin(the)*sin(phi)+cos(psi)*cos(phi) sin(psi)*sin(the)*cos(phi)-cos(psi)*sin(phi);
           -sin(the)         cos(the)*sin(phi)                            cos(the)*cos(phi)];

    % Support for X-configuration
    if (~evalin('base','quadModel.plusConfig'))
        Rz = [ sqrt(2)/2, sqrt(2)/2, 0;
               -sqrt(2)/2,sqrt(2)/2, 0;
                       0,          0, 1];
        Nr=Rz*Nr;
        Er=Rz*Er;
        Wr=Rz*Wr;
        Sr=Rz*Sr;
        mN = Rz*mN;
        mE = Rz*mE;
        mW = Rz*mW;
        mS = Rz*mS;
        bNS = Rz*bNS;
        bEW = Rz*bEW;
        Top = Rz*Top;
        Bot = Rz*Bot;
        NEB = Rz*NEB;
        NWB = Rz*NWB;
        SWB = Rz*SWB;
        SEB = Rz*SEB;
    end

    % Bring in and rotate body frame velocity vector
    U = A(:,7);
    V = A(:,8);
    W = A(:,9);
    Vi = zeros(length(A),3);
    MvMax= max(sqrt(U.^2+V.^2+W.^2));
    Vb = 3/MvMax*[U, V, W].'; % Scale velocity vector points for plotting
    Vi(1,:) = R*Vb(:,1); % Rotate velocity vector to inertial frame for plotting
    % Angular Velocity Vector
    P = A(:,1);
    Q = A(:,2);
    Rw = A(:,3);
    omi = zeros(length(A),3); % Initialize omega inertiaF points array (L(A)x3)
    MombMax = max(sqrt(P.^2+Q.^2+Rw.^2)); % Calculate max magnitude of omb
    omb = 3/MombMax*[P,Q,Rw].'; % Store and scale current omega bodyF
    omi(1,:) = R*omb(:,1); % Rotate omegab to inertiaF store in omegai array

    phi = A(:,4);
    the= A(:,5);
    psi = A(:,6);

    % Next we run through the points in the vector A as an
    % animation until we get a command to stop
    minX = min(A(:,10))*3.28; % meters to feet
    minY = min(A(:,11))*3.28;
    %minZ = min(A(:,12))*3.28;
    maxX = max(A(:,10))*3.28;
    maxY = max(A(:,11))*3.28;
    maxZ = max(A(:,12))*3.28;
    colors = jet(length(A(:,10))); % color the path for time info

    % Precompute R
    % ROTATION MATRIX BELOW --- ZYX ROTATION
    R = cell(length(A),1); % create empty cell array for R matrices
    for i = 1:length(A)
    R{i,1} = [cos(psi(i))*cos(the(i)) cos(psi(i))*sin(the(i))*sin(phi(i))-sin(psi(i))*cos(phi(i)) cos(psi(i))*sin(the(i))*cos(phi(i))+sin(psi(i))*sin(phi(i));
              sin(psi(i))*cos(the(i)) sin(psi(i))*sin(the(i))*sin(phi(i))+cos(psi(i))*cos(phi(i)) sin(psi(i))*sin(the(i))*cos(phi(i))-cos(psi(i))*sin(phi(i));
              -sin(the(i))         cos(the(i))*sin(phi(i))                            cos(the(i))*cos(phi(i))];
    end
    % set up a j vector to straighten out skipping issues
%     if (skipFlag == 1)
%         jSkip = 1:frameSkipVal:j;
%     end
    while ( strcmp(get(hObject,'String'),'Stop'))
        handles = guidata(gcbo);
        j = handles.j;

        guidata(hObject,handles); % Update handles data
        Vi(j,:) = R{j,1}*Vb(:,j);

        NrR = R{j,1}*Nr;
        ErR = R{j,1}*Er;
        WrR = R{j,1}*Wr;
        SrR = R{j,1}*Sr;
        bNSR = R{j,1}*bNS;
        bEWR = R{j,1}*bEW;
        TopR = R{j,1}*Top;
        BotR = R{j,1}*Bot;
        % Rotate body parts Via Initialized R
        mNr = R{j,1}*mN;
        mEr = R{j,1}*mE;
        mWr = R{j,1}*mW;
        mSr = R{j,1}*mS;
        NEBR = R{j,1}*NEB;
        NWBR = R{j,1}*NWB;
        SWBR = R{j,1}*SWB;
        SEBR = R{j,1}*SEB;

        % Plot the quad rotation and ang. velocity and inertial frame velocity vector
        axes(handles.axes1)
        plot3(bNSR(1,:),bNSR(2,:),bNSR(3,:),'b','LineWidth',3)
        hold on
        plot3(bEWR(1,:),bEWR(2,:),bEWR(3,:),'b','LineWidth',3)
        plot3(NrR(1,:),NrR(2,:),NrR(3,:),'g')
        plot3(ErR(1,:),ErR(2,:),ErR(3,:),'k')
        plot3(WrR(1,:),WrR(2,:),WrR(3,:),'k')
        plot3(SrR(1,:),SrR(2,:),SrR(3,:),'g')
        plot3(mNr(1,:),mNr(2,:),mNr(3,:),'k','LineWidth',4) % Motor
        plot3(mEr(1,:),mEr(2,:),mEr(3,:),'k','LineWidth',4) % Motor
        plot3(mWr(1,:),mWr(2,:),mWr(3,:),'k','LineWidth',4) % Motor
        plot3(mSr(1,:),mSr(2,:),mSr(3,:),'k','LineWidth',4) % Motor
        fill3(TopR(1,:),TopR(2,:),TopR(3,:),'r'); alpha(0.8);
        fill3(BotR(1,:),BotR(2,:),BotR(3,:),'g'); alpha(0.8);
        grey = [0.5 0.5 0.5];
        ne  = fill3(NEBR(1,:),NEBR(2,:),NEBR(3,:),'c'); alpha(ne,0.8); % North East surface
        nw  = fill3(NWBR(1,:),NWBR(2,:),NWBR(3,:),grey); alpha(nw,0.8); % North West surface
        sw  = fill3(SWBR(1,:),SWBR(2,:),SWBR(3,:),grey); alpha(sw,0.8); % South West surface
        se  = fill3(SEBR(1,:),SEBR(2,:),SEBR(3,:),grey); alpha(se,0.8); % South East surface
        xlabel('X')
        ylabel('Y')
        zlabel('Z')
        xlim([-3 3])
        ylim([-3 3])
        zlim([-3 3])
        % set(gca,'Color',[215/255 244/255 247/255]) % This is a technique
        % that can be used for setting the plot background color
        handles.AZval = get(handles.AZslider,'Value')*360;
        handles.ELval = get(handles.ELslider,'Value')*90;
        view(handles.axes1,handles.AZval,handles.ELval)
%         P = A(j,1); Q = A(j,2); Rw = A(j,3);
%         M = sqrt(P^2+Q^2+Rw^2);
%         omb = 3/M*[P,Q,Rw].';%/M; % Scaling of angular velocity vector
        omi(j,:) = R{j,1}*omb(:,j);
        qp1 = quiver3(0,0,0,omi(j,1),omi(j,2),omi(j,3),'r');
        qp2 = quiver3(0,0,0,Vi(j,1),Vi(j,2),Vi(j,3),'k');
        axis square
        grid on
        hold off
        drawnow

        % Plot the three dimensional trajectory of the box
        if (j==1)
            cla(handles.axes2)
        else if (handles.skipFlag==1)
                cla(handles.axes2)
                X = A(1:frameSkipVal:j,10)*3.28;
                Y = A(1:frameSkipVal:j,11)*3.28;
                Z = A(1:frameSkipVal:j,12)*3.28;
                axes(handles.axes2)
                hold on
                scatter3(X,Y,Z,36,colors(1:frameSkipVal:j,:));
            end
        end
        % Ordinary plotting sequence (one frame at a time using hold on to
        % create persistance
        X = A(j,10)*3.28; Y = A(j,11)*3.28; Z = A(j,12)*3.28;
        axes(handles.axes2)
        hold on
        scatter3(X,Y,Z,36,colors(j,:));
        if (j == 1 || handles.skipFlag==1)
            fill3([minX-1 maxX+1 maxX+1 minX-1],...
                  [minY-1 minY-1 maxY+1 maxY+1],...
                  [0 0 0 0],'g'); % make a plane to represent the ground (Z = 0)
            alpha(0.7); % Makes the ground "see-through"
            handles.skipFlag = 0;
            guidata(hObject,handles)
        end
        xlabel('X (ft)')
        ylabel('Y (ft)')
        zlabel('Z (ft)')
        xlim([minX-1 maxX+1])
        ylim([minY-1 maxY+1])
        zlim([-0.1 maxZ+1]) % Keep ground within view
        % These calls are necessarry to avoid stutter in the frames if adjusting slider while animation is running
        handles.AZval = get(handles.AZslider,'Value')*360;
        handles.ELval = get(handles.ELslider,'Value')*90;
        view(handles.axes2,handles.AZval,handles.ELval) % Set the current vie
        axis square
%         axis equal
%         xlimz = xlim;
%         ylimz = ylim;
%                     fill3([xlimz(1) xlimz(2) xlimz(2) xlimz(1)],...
%                   [ylimz(1) ylimz(1) ylimz(2) ylimz(2)],...
%                   [0 0 0 0],'k'); % make a plane to represent the ground (Z = 0)
%             alpha(0.9); % Makes the ground "see-through"
        grid on
        drawnow
        hold off
%         if (firstRunFlag==1)
%             firstRunFlag = 0;
%         end
        if (j == size(A,1)) % If we've reached the last group of data
            handles.j = 1; % Reset j to 1
            set(hObject,'String','Start'); % Reset button state
            guidata(hObject,handles);
            return % exit the function

        else if (j+frameSkipVal<size(A,1)) % If the next index won't exceed total frames...
        j = j+frameSkipVal; % increase the index by frame skip value
            else
            j = size(A,1); % We're almost out of data so to the last value for the final plot
            end
        end
        set(handles.simTime,'String',num2str(tout(j))); % update time display
        handles = guidata(gcbo); % update handles structure again...
        if (handles.skipFlag == 1)
            guidata(hObject,handles) % this nightmare is used to catch skip events which otherwise tend to go unheeded.
            return
        else
        handles.j = j; % store j as a handles value for persistance
        guidata(hObject,handles);
        end
    end
    else if strcmp(get(hObject,'String'),'Stop')
        set(hObject,'String','Start');
        end
end

% --- Executes during object creation, after setting all properties.
function simTime_CreateFcn(hObject, eventdata, handles)
% hObject    handle to simTime (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% --- Executes during object creation, after setting all properties.
function axes2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to axes2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate axes2

% --- Executes on button press in XYview.
function XYview_Callback(hObject, eventdata, handles)
% hObject    handle to XYview (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
set(handles.ELslider,'Value',1);
set(handles.AZslider,'Value',0);
view(handles.axes1,0,90);
view(handles.axes2,0,90);
drawnow
guidata(hObject,handles);

% --- Executes on button press in XZview.
function XZview_Callback(hObject, eventdata, handles)
% hObject    handle to XZview (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
set(handles.ELslider,'Value',0);
set(handles.AZslider,'Value',0);
view(handles.axes1,0,0);
view(handles.axes2,0,0);
drawnow
guidata(hObject,handles);

% --- Executes on button press in YZview.
function YZview_Callback(hObject, eventdata, handles)
% hObject    handle to YZview (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
set(handles.ELslider,'Value',0);
set(handles.AZslider,'Value',90/360);
view(handles.axes1,90,0);
view(handles.axes2,90,0);
drawnow
guidata(hObject,handles);

% --- Executes on button press in defaultView.
function defaultView_Callback(hObject, eventdata, handles)
% hObject    handle to defaultView (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
set(handles.ELslider,'Value',25/90);
set(handles.AZslider,'Value',35/360)
view(handles.axes1,35,25);
view(handles.axes2,35,25);
drawnow
guidata(hObject,handles);

% --- Executes on button press in timeSkipButton.
function timeSkipButton_Callback(hObject, eventdata, handles)
% hObject    handle to timeSkipButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
set(handles.startStop,'String','Start');
tout = evalin('base', 'tout');
TimeSkip = str2double(get(handles.timeSkipEditBox,'String'));
handles.j = find(tout<=TimeSkip, 1,'last');
handles.skipFlag = 1;
% if handles.doubleRunFlag == 0
%     timeSkipButton_Callback(hObject, eventdata, handles) % run it again
% end
% handles.doubleRunFlag = 1;
guidata(hObject,handles);
% yes i'm seriously resorting to this...

function timeSkipEditBox_Callback(hObject, eventdata, handles)
% hObject    handle to timeSkipEditBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of timeSkipEditBox as text
%        str2double(get(hObject,'String')) returns contents of timeSkipEditBox as a double
handles.TimeSkip = str2double(get(hObject,'String'));
handles.skipFlag = 0;
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function timeSkipEditBox_CreateFcn(hObject, eventdata, handles)
% hObject    handle to timeSkipEditBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
handles.TimeSkip = 0;
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function figure1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

repo: Quad-Sim

path: github-repos/MATLAB/dch33-Quad-Sim/Quad-Sim-961bc69d4939c8d0661eeb9820de668994262f65/Quadcopter Dynamic Modeling and Simulation/Graphical User Interfaces/QuadAnim4.m


text:
function ConvHDRvtoLDRv(hdrv, filenameOutput, fstops, ldrv_gamma, ldrv_quality, ldrv_video_profile)
%
%
%        ConvHDRvtoLDRv(hdrv, filenameOutput, fstops, ldrv_gamma, ldrv_quality, ldrv_video_profile)
%
%
%
%        Input:
%           -hdrv:
%           -filenameOutput:
%           -fstops:
%           -ftops:
%           -ldrv_gamma:
%           -ldrv_quality:
%           -ldrv_video_profile:
%
%
%     Copyright (C) 2016  Francesco Banterle
%
%     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(~exist('ldrv_gamma', 'var'))
    ldrv_gamma = 2.2;
end

if(~exist('ldrv_quality', 'var'))
    ldrv_quality = 95;
end

if(~exist('ldrv_video_profile', 'var'))
    ldrv_video_profile = 'MPEG-4';
end

if(ldrv_gamma < 0)
    bsRGB = 1;
else
    bsRGB = 0;
end

if(isempty(fstops))
    errir('ConvHDRvtoLDRv: fstops cannot be empty!');
end

name = RemoveExt(filenameOutput);
ext = fileExtension(filenameOutput);

bVideo = 0;
writerObj = 0;

if(strcmp(ext, 'avi') == 1 | strcmp(ext, 'mp4') == 1)
    bVideo = 1;
    writerObj = VideoWriter(filenameOutput, ldrv_video_profile);
    writerObj.FrameRate = hdrv.FrameRate;
    writerObj.Quality = ldrv_quality;
    open(writerObj);
end

if(bVideo == 0)
    mkdir([name,'_img']);
end

hdrv = hdrvopen(hdrv, 'r');

n = length(fstops);

for i=1:hdrv.totalFrames
    disp(['Processing frame ', num2str(i)]);
    [frame, hdrv] = hdrvGetFrame(hdrv, i);

    frame(frame < 0) = 0;

    j = mod(i, n) + 1;

    frameOut = frame * 2^fstops(j);

    %Gamma/sRGB encoding
    if(bsRGB)
        frameOut = ClampImg(ConvertRGBtosRGB(frameOut, 0), 0, 1);
    else
        frameOut = ClampImg(GammaTMO(frameOut, ldrv_gamma, 0, 0), 0, 1);
    end

    %Storing
    if(bVideo)
        writeVideo(writerObj, frameOut);
    else
        nameOut = [name, '_img/frame_', sprintf('%.10d',i), '.', ext];
        imwrite(frameOut, nameOut);

        nameOut = [name, '_img/frame_', sprintf('%.10d',i), '.exr'];
        hdrimwrite(frame, nameOut);
    end
end

hdrvclose(hdrv);

end

author: banterle


text:
function metScores = compareAllMets(Cmodel, Tmodel)
% Compares and scores similarity between all pairwise
% combinations of metabolites between two models. It is not however used
% by the main workflow.
%
% USAGE:
%
%    metScores = compareAllMets(Cmodel,Tmodel
%
% INPUTS:
%    Cmodel:      Comparison model.
%    Tmodel:      Template model.
%
% OUTPUTS:
%    metScore:    Matrix of scores of all pairwise comparisons.
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

nCmets = length(Cmodel.mets) ; % Declare variables and scoring structure.
nTmets = length(Tmodel.mets) ;

% Score values.
ScoreVal = struct('ID', [10,-1], ...
                  'compartment', [0,-100], ...
                  'name', [10,-1], ...
                  'formula', [10,-5], ...
                  'charge', [1,-1], ...
                  'KEGGID', [10,-5], ...
                  'SEEDID', [10,-5], ...
                  'ChEBIID', [10,-5], ...
                  'PubChemID', [10,-5], ...
                  'InChIString', [10,-5]) ;

% Score matrix.
metScores = zeros(nCmets, nTmets) ;

% Compartments of every met in T, used later for comparison.
tComp = cell(nTmets, 1) ;
for iMet = 1:nTmets
    tComp{iMet} = Tmodel.mets{iMet}(end - 2:end) ;
end

%% Compare cMet to tMets
h = waitbar(0, 'Comparing metabolites') ;
for iMet = 1:nCmets
    % metsID Match. There is always one ID.
    % Match name, removing '[.]' and add regular expression info.
    if ~isempty(regexp(Cmodel.mets{iMet}, '\[.\]', 'once'))
        name =  strcat('(\||\<)', Cmodel.mets{iMet}(1:end - 3), '(\>|\[|_|\|)') ;
    else
        name = strcat('(\||\<)', Cmodel.mets{iMet}, '(\>|\[|_|\|)') ;
    end
    % Match against IDs for all models.
    match = regexpi(Tmodel.metID, name) ;
    match = ~cellfun(@isempty, match) ;
    % Give scores.
    metScores(iMet, match) = metScores(iMet, match) + ScoreVal.ID(1) ;
    metScores(iMet, ~match) = metScores(iMet, ~match) + ScoreVal.ID(2) ;

    % Match against compartment for all models.
    match = strcmp(tComp, Cmodel.mets{iMet}(end - 2:end)) ;
    % Give scores.
    metScores(iMet, match) = metScores(iMet, match) + ScoreVal.compartment(1) ;
    metScores(iMet, ~match) = metScores(iMet, ~match) + ScoreVal.compartment(2) ;

    % metNames Match. There can be multiple names.
    fullname = Cmodel.metNames{iMet} ;
    pipePos = [0 strfind(fullname, '|') length(fullname) + 1] ;
    matchSum = zeros(length(Tmodel.mets), 1) ;
    for i = 1:length(pipePos) - 1
        name = fullname(pipePos(i) + 1:pipePos(i + 1) - 1) ;
        name = strcat('(\||\<)', name, '(\>|\|)') ;
        match = regexpi(Tmodel.metNames, name) ;
        match = ~cellfun('isempty', match) ;
        matchSum = matchSum + match ;
    end
    match = logical(matchSum) ;
    % Give scores.
    metScores(iMet, match) = metScores(iMet, match) + ScoreVal.name(1) ;
    metScores(iMet, ~match) = metScores(iMet, ~match) + ScoreVal.name(2) ;

    % metFormula Match. There will always be one Formula.
    name = strcat('(\||\<)', Cmodel.metFormulas{iMet}, '(\>|\|)') ;
    match = regexpi(Tmodel.metFormulas, name) ;
    match = ~cellfun(@isempty, match) ;
    % Give scores.
    metScores(iMet, match) = metScores(iMet, match) + ScoreVal.formula(1) ;
    metScores(iMet, ~match) = metScores(iMet, ~match) + ScoreVal.formula(2) ;

    % metCharge match. There will always be a charge.
    for i = 1:length(Tmodel.mets)
        if Cmodel.metCharge(iMet) == Tmodel.metCharge(i)
            metScores(iMet, i) = metScores(iMet, i) + ScoreVal.charge(1) ;
        else
            metScores(iMet, i) = metScores(iMet, i) + ScoreVal.charge(2) ;
        end
    end

    % metKEGGID Match. There may multiple KEGGIDs or no KEGGID.
    if ~isempty(Cmodel.metKEGGID{iMet})
        fullname = Cmodel.metKEGGID{iMet} ;
        pipePos = [ 0 strfind(fullname, '|') length(fullname) + 1 ] ;
        matchSum = zeros(length(Tmodel.mets), 1) ;
        for i = 1:length(pipePos) - 1
            name = fullname(pipePos(i) + 1:pipePos(i + 1) - 1) ;
            match = regexpi(Tmodel.metKEGGID, name) ;
            match = ~cellfun('isempty', match) ;
            matchSum = matchSum + match ;
        end
        match = logical(matchSum) ;
        % Give scores.
        metScores(iMet, match) = metScores(iMet, match) + ScoreVal.KEGGID(1) ;
        metScores(iMet, ~match) = metScores(iMet, ~match) + ScoreVal.KEGGID(2) ;
    end

    % metSEEDID Match. There may multiple SEEDIDs or no SEEDID.
    if ~isempty(Cmodel.metSEEDID{iMet})
        fullname = Cmodel.metSEEDID{iMet} ;
        pipePos = [ 0 strfind(fullname, '|') length(fullname) + 1] ;
        matchSum = zeros(length(Tmodel.mets), 1) ;
        for i = 1:length(pipePos) - 1
            name = fullname(pipePos(i) + 1:pipePos(i +1 ) - 1) ;
            match = regexpi(Tmodel.metSEEDID,name) ;
            match = ~cellfun('isempty', match) ;
            matchSum = matchSum + match ;
        end
        match = logical(matchSum) ;
        % Give scores.
        metScores(iMet, match) = metScores(iMet, match) + ScoreVal.SEEDID(1) ;
        metScores(iMet, ~match) = metScores(iMet, ~match) + ScoreVal.SEEDID(2) ;
    end
    % Update waitbar.
    waitbar(iMet / nCmets, h)
end
close(h)

%% Normalize scores.
% remove negative scores
metScores(metScores < 0) = 0 ;

if max(max(metScores)) > 0 % Avoid NaN.
    metScores = metScores - min(min(metScores)) ;
    metScores = metScores ./ max(max(metScores)) ;
else
    metScores = zeroes(size(metScores)) ;
end

% Find best matches.
[bestMatch, bestMatchIndex] = max(metScores, [], 2) ;

%% Graph data.
hold on
figure(1)
subplot(2,1,1)
hist(metScores(:), 100) % Converting to column matrix is faster.
title('Normalized score frequency')
xlabel('Score')
ylabel('Frequency')
subplot(2, 2, 3)
hist(bestMatch, 100)
title('Frequency of scores of best match per metabolite')
xlabel('Score')
ylabel('Frequency');
subplot(2, 2, 4)
imagesc(metScores)
colormap('bone')
title('Matching scores')
xlabel('Template met')
ylabel('Compared met')
hold off

author: opencobra

repo: cobratoolbox


text:
function b = Theil_Sen_Regress(x,y)

[N c]=size(x);

Comb = combnk(1:N,2);
deltay=diff(y(Comb),1,2);
deltax=diff(x(Comb),1,2);

theil=diff(y(Comb),1,2)./diff(x(Comb),1,2);
b=median(theil);

%% Example code below !!

% % %  x = (-15:25)';
% % % y =2*x + normrnd(0,3,41,1);
% % % y([38 40 41]) = 0;
% % % y([1 3])=20;
% % %
% % % bls = regress(y,[ones(41,1) x]);
% % % [N c]=size(x);
% % % Comb = combnk(1:N,2);
% % % deltay=diff(y(Comb),1,2);
% % % deltax=diff(x(Comb),1,2);
% % % theil=diff(y(Comb),1,2)./diff(x(Comb),1,2);
% % % b=median(theil);
% % %
% % % scatter(x,y,'filled'); grid on; hold on
% % % plot(x,bls(1)+bls(2)*x,'r','LineWidth',2);
% % % plot(x,b*x,'g','LineWidth',3);
% % % plot(x,2*x,'k','LineWidth',2);
% % % legend('Data','Ordinary Least Squares','Theil-Senn Estimator','Real slope')


text:
%TR2RPY	Convert a homogeneous transform matrix to roll/pitch/yaw angles
%
%	[A B C] = TR2RPY(TR) returns a vector of Euler angles
%	corresponding to the rotational part of the homogeneous transform TR.
%
%	See also  RPY2TR, TR2EUL

%	Copright (C) Peter Corke 1993
function rpy = tr2rpy(m)

	rpy = zeros(1,3);

	if abs(m(1,1)) < eps & abs(m(2,1)) < eps,
		rpy(1) = 0;
		rpy(2) = atan2(-m(3,1), m(1,1));
		rpy(3) = atan2(-m(2,3), m(2,2));
	else,
		rpy(1) = atan2(m(2,1), m(1,1));
		sp = sin(rpy(1));
		cp = cos(rpy(1));
		rpy(2) = atan2(-m(3,1), cp * m(1,1) + sp * m(2,1));
		rpy(3) = atan2(sp * m(1,3) - cp * m(2,3), cp*m(2,2) - sp*m(1,2));
	end

repo: ekfmonoslam


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

function ip2=preprocess(ip1,fs1,intval1)
[rw,cm]=size(ip1);
if rw==1
    ip1=ip1';
else ip1=ip1;
end;
ip1=ip1(:,1);                %to make one dimentional
len1=length(ip1);
i = 1;
while abs(ip1(i)) < 0.002 && i<48000% Silence detection
i = i + 1;
end
ip1(1:i) = [];
ip1=ip1-mean(ip1);
ip1=ip1*2;                 % amplification
len1=length(ip1);
i=fs1*intval1;
ip2=zeros(i,1);        % to make uniform dimention size
if len1<i
    ip2(1:len1)=ip1(1:len1);
else
    ip2(1:i)=ip1(1:i);
end;


text:
function range= visutil_commonRangeForGA(erp, varargin)
%VISUTIL_COMMONRANGEFORGA - Determine plotting range for grand average ERPs
%
%Synposis:
% RANGE= visutil_commonRangeForGA(ERP, <OPT>)
%
%Input:
% ERP: cell of size [N x C] containing ERPs of N subjects and C conditions,
%      plotting ranges are determined for each condition separately
% OPT: struct or property/value list of optional properties:
%  .CLabERP, CLabScalp - channels to be considered for ERP/scalp plots,
%                        default: '*' (all channels)
%  .IvalERP, IvalScalp - interval for ERP/scalp plots to be considered,
%                        default: [] (whole ERP)
%  .SymERP, SymScalp   - make ERP/scalp range symmetric,
%                        default: false/true
%
%Output:
% RANGE: struct containing the two fields:
%  .erp   - plotting range for single channel ERPs
%  .scalp - plotting range for scalp maps

props= {'CLabERP'         '*'    'CELL{CHAR}|CHAR';
        'CLabScalp'       '*'    'CELL{CHAR}|CHAR';
        'IvalERP'         []     'DOUBLE[2]';
        'IvalScalp'       []     'DOUBLE';
        'SymERP'          0      'BOOL';
        'SymScalp'        1      'BOOL'
        'NiceRangeERP'    0      'DOUBLE';
        'EnlageRangeERP'  0.02   'DOUBLE';
       };

if nargin==0,
  range= props; return
end

misc_checkType(erp, 'CELL{STRUCT(x fs clab)}');
opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

if opt.NiceRangeERP && isdefault.EnlageRangeERP,
  opt.EnlageRangeERP= 0;
end

range.erp= [inf -inf];
if ~isempty(opt.IvalScalp),
  range.scalp= [inf -inf];
end

for ff= 1:size(erp,2),
  erp_ga= proc_grandAverage(erp{:,ff});

  tmp= proc_selectChannels(erp_ga, opt.CLabERP);
  if ~isempty(opt.IvalERP),
    tmp= proc_selectIval(tmp, opt.IvalERP);
  end
  range.erp(1)= min(range.erp(1), min(tmp.x(:)));
  range.erp(2)= max(range.erp(2), max(tmp.x(:)));

  if ~isempty(opt.IvalScalp),
    tmp= proc_jumpingMeans(erp_ga, opt.IvalScalp);
    tmp= proc_selectChannels(tmp, opt.CLabScalp);
    range.scalp(1)= min(range.scalp(1), min(tmp.x(:)));
    range.scalp(2)= max(range.scalp(2), max(tmp.x(:)));
  end

end

if opt.SymERP,
  range.erp= [-1 1]*max(abs(range.erp));
end
if opt.EnlageRangeERP>0,
  range.erp= range.erp + [-1 1]*opt.EnlageRangeERP*diff(range.erp);
end
if opt.NiceRangeERP>0,
  res= opt.NiceRangeERP;
  range.erp(1)= floor(res*range.erp(1))/res;
  range.erp(2)= ceil(res*range.erp(2))/res;
end
if ~isempty(opt.IvalScalp) && opt.SymScalp,
  range.scalp= [-1 1]*max(abs(range.scalp));
end

author: bbci


text:
%SingleCell
%
% Author: Derek Macklin, macklin@stanford.edu
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 7/6/2011
classdef SingleCell
    methods (Static = true)
        function run(simBatchDir, simNum, fileName, methodsToRun)
            import edu.stanford.covert.cell.sim.analysis.SingleCell;
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;
            import edu.stanford.covert.cell.sim.util.PlotUtil;

            %options
            if nargin < 1 || isempty(simBatchDir)
                simBatchDir = SimulationDiskUtil.getLatestSimulationGroup();
            end

            if nargin < 2 || isempty(simNum)
                simNum = 1;
            end

            %simulation
            [~, ~, sim] = SimulationDiskUtil.getSimulation([simBatchDir filesep num2str(simNum)]);

            %evaluate each method and optionally save at the result as
            %pdf files
            metaData = meta.class.fromName('edu.stanford.covert.cell.sim.analysis.SingleCell');
            for i = 1:numel(metaData.Methods)
                %only plot methods matching the signature
                %  plotFcn(sim, simBatchDir, simNum, figHandle)
                if ...
                        numel(metaData.Methods{i}.InputNames) < 4 || ~isequal(metaData.Methods{i}.InputNames(1:4), {'sim'; 'simBatchDir'; 'simNum'; 'figHandle'}) || ...
                        numel(metaData.Methods{i}.OutputNames) < 1 || ~isequal(metaData.Methods{i}.OutputNames(1), {'figData'})
                    continue;
                end

                %get method name
                methodName = metaData.Methods{i}.Name;
                MethodName = [upper(methodName(1)) methodName(2:end)];
                if nargin >= 4 && ~ismember(metaData.Methods{i}.Name, methodsToRun)
                    continue;
                end

                %run plot method
                [~, figHandle] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle();
                try
                    argout = cell(size(metaData.Methods{i}.OutputNames));
                    [argout{:}] = SingleCell.(methodName)(sim, simBatchDir, simNum, figHandle);

                    % save plot
                    if nargin >= 3
                        figData = argout{1}; %#ok<NASGU>
                        save([fileName '-' MethodName '.mat'], '-struct', 'figData');

                        orient(figHandle, 'portrait');
                        saveas(figHandle, [fileName '-' MethodName '.pdf']);
                    end

                    % save table
                    if numel(argout) > 1
                        if nargin >= 3
                            PrintUtil.printToFile(argout{2}, argout{3}, [fileName '-' MethodName '.xls'], MethodName);
                        else
                            PrintUtil.printToStdIO(argout{2}, argout{3});
                        end
                    end
                catch exception
                    warning('WholeCell:warning', 'Unable to make %s plot:\n%s', methodName, exception.getReport());
                end

                if nargin >= 2
                    close(figHandle);
                    clear figHandle;
                end
                clear argout figData;
            end
        end

        function figData = growth(sim, simBatchDir, simNum, figHandle)
            import edu.stanford.covert.cell.sim.analysis.ChromosomeSpaceTimePlot;
            import edu.stanford.covert.cell.sim.analysis.SingleCell;
            import edu.stanford.covert.cell.sim.util.PlotUtil;

            %% get data
            [time, mass, rnas, proteins, ...
                dntps, ntps, aminoAcids, ...
                growth, reactionFlux, ...
                enzymeComplex, enzymeMonomer, enzymeRNA, enzymeGene, ...
                ploidy, superhelicalDensity] = ...
                SingleCell.calcGrowthData(sim, simBatchDir, simNum);

            %% layout figure, plot data
            clf(figHandle);

            options = struct();
            options.titleStr = {'Composition'; 'Growth'; 'Organization'};
            options.position = [0.12 0.52 0.87 0.4];
            options.xdata = time;
            options.ydata = {{
                permute(mass, [1 3 2 4])
                permute(rnas, [1 3 2 4])
                permute(proteins, [1 3 2 4])
                permute(dntps, [1 3 2 4])
                permute(ntps, [1 3 2 4])
                permute(aminoAcids, [1 3 2 4])
                }
                {
                permute(growth, [1 3 2 4])
                permute(reactionFlux, [1 3 2 4])
                permute(enzymeComplex, [1 3 2 4])
                permute(enzymeMonomer, [1 3 2 4])
                permute(enzymeRNA, [1 3 2 4])
                permute(enzymeGene, [1 3 2 4])
                }
                {
                permute(ploidy, [1 3 2 4])
                permute(superhelicalDensity, [1 3 2 4])
                []
                }};
            options.ylabelStr = {{
                'Mass (fg)'
                'RNAs'
                'Proteins'
                'dNTPs'
                'NTPs'
                {'Amino' 'acids'}
                }
                {
                {'Growth' '(fg h^{-1})'}
                {sprintf('%s', rxnId) 'flux' '(MHz)'}
                {sprintf('%s', rxnId) 'complex'}
                {sprintf('%s', rxnId) 'monomer'}
                {sprintf('%s', rxnId) 'mRNA'}
                sprintf('{\\it %s}', rxnId)
                }
                {
                {'Copy', 'number'}
                {'Super-' 'helicity'}
                {'Polymerase' 'position'}
                }
                };
            axesHandles = PlotUtil.multiElementPlot(figHandle, {3 * ones(6, 1); 3 * ones(6, 1); [3 3 4*3+(4-1)*1]}, time([1 end])', ...
                options);
            ChromosomeSpaceTimePlot.plotSpaceTime(axesHandles{3}(3), simBatchDir, simNum, [], false);

            clear options;

            %% Format Y-axis
            PlotUtil.labelSubplots([axesHandles{1}(1); axesHandles{2}(1); axesHandles{3}(1)]);

            %% collect data
            figData = struct;
            figData.time = time;
            figData.mass = mass;
            figData.rnas = rnas;
            figData.proteins = proteins;
            figData.dntps = dntps;
            figData.ntps = ntps;
            figData.aminoAcids = aminoAcids;
            figData.growth = growth;
            figData.reactionFlux = reactionFlux;
            figData.enzymeComplex = enzymeComplex;
            figData.enzymeMonomer = enzymeMonomer;
            figData.enzymeRNA = enzymeRNA;
            figData.enzymeGene = enzymeGene;
            figData.ploidy = ploidy;
            figData.superhelicalDensity = superhelicalDensity;
        end

        function [time, mass, rnas, proteins, ...
                dntps, ntps, aminoAcids, ...
                growth, reactionFlux, ...
                enzymeComplex, enzymeMonomer, enzymeRNA, enzymeGene, ...
                ploidy, superhelicalDensity] = ...
                calcGrowthData(sim, simBatchDir, simNum, rxnId)

            import edu.stanford.covert.cell.sim.util.MacromoleculeUtil;
            import edu.stanford.covert.cell.sim.util.SimulationEnsemble;
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;
            import edu.stanford.covert.cell.sim.analysis.CellOverview;
            import edu.stanford.covert.util.ConstantUtil;

            comp = sim.compartment;
            g = sim.gene;

            c = sim.state('Chromosome');
            m = sim.state('Metabolite');
            rna = sim.state('Rna');
            pm = sim.state('ProteinMonomer');
            pc = sim.state('ProteinComplex');

            sc = sim.process('DNASupercoiling');
            met = sim.process('Metabolism');

            %% get data
            stateNames = {
                'mass'              'Mass'
                'ploidy'            {'Chromosome' 'Copy Number'}
                'rnas'              'RNA'
                'proteins'          'Proteins'
                'lipids'            'Lipids'
                'growth_rate'       'Growth'
                'amino_acids'       {'Amino' 'Acids'}
                'ntps'              'NTPs'
                };
            ensemble = SimulationEnsemble(simBatchDir, stateNames, [], simNum);

            time = ensemble.stateData.time / 3600;
            mass = ensemble.stateData.values(ensemble.getPropertyIndices('mass'), :, :, :) * sim.state('Mass').cellInitialDryWeight / (1-sim.state('Mass').fractionWetWeight) * 1e15;
            ploidy = ensemble.stateData.values(ensemble.getPropertyIndices('ploidy'), :, :, :);
            rnas = ensemble.stateData.values(ensemble.getPropertyIndices('rnas'), :, :, :);
            proteins = ensemble.stateData.values(ensemble.getPropertyIndices('proteins'), :, :, :);
            %lipids = ensemble.stateData.values(ensemble.getPropertyIndices('lipids'), :, :, :);
            growth = ensemble.stateData.values(ensemble.getPropertyIndices('growth_rate'), :, :, :) * 3600 * sim.state('Mass').cellInitialDryWeight / (1-sim.state('Mass').fractionWetWeight)  * 1e15;

            clear ensemble;

            if nargin < 4
                rxnId = 'AckA';
            end
            rxnIdx = met.reactionIndexs(rxnId);
            enzId = met.enzymeWholeCellModelIDs{met.reactionCatalysisMatrix(rxnIdx, :) ~= 0};
            [~, ~, ~, ~, ~, ...
                geneIdxs, ~, ~, ~, ...
                ~, ~, ~, ~, ...
                matureRnaIdxs, ~, ~, ~, ...
                monomerIdxs, ~, ~, ~, ...
                complexIdxs, ~, ~, ~] = ...
                MacromoleculeUtil.getMacromoleculeIndexsIDsNames(enzId, sim);
            stateNames = {
                'Metabolite'        'counts'              ':'  ':'
                'MetabolicReaction' 'fluxs'               rxnIdx ':'
                'Chromosome'        'polymerizedRegions'  ':'  ':'
                'Chromosome'        'linkingNumbers'       ':'  ':'
                'Rna'               'counts' rna.matureIndexs(matureRnaIdxs)         comp.cytosolIndexs
                'ProteinMonomer'    'counts' [pm.matureIndexs([monomerIdxs; sc.enzymeGlobalIndexs(sc.enzymeIndexs_topoI)]); pm.boundIndexs([monomerIdxs; sc.enzymeGlobalIndexs(sc.enzymeIndexs_topoI)])] comp.cytosolIndexs
                'ProteinComplex'    'counts' [pc.matureIndexs([complexIdxs; sc.enzymeGlobalIndexs([sc.enzymeIndexs_gyrase; sc.enzymeIndexs_topoIV])]); pc.boundIndexs([complexIdxs; sc.enzymeGlobalIndexs([sc.enzymeIndexs_gyrase; sc.enzymeIndexs_topoIV])])] comp.cytosolIndexs
                };
            states = SimulationEnsemble.load(simBatchDir, stateNames, [], [], 1, 'extract', simNum);
            tmp = full(states.Metabolite.counts);
            %smallMolecules = sum(sum(tmp(setdiff(1:end, m.lipidIndexs), [comp.cytosolIndexs; comp.membraneIndexs], :), 2), 1);
            dntps = cat(3, NaN(4, 1), tmp(m.dntpIndexs, 1, :));
            ntps = cat(3, NaN(4, 1), tmp(m.ntpIndexs, 1, :));
            %nmps = tmp(m.nmpIndexs(1:20), comp.cytosolIndexs, :);
            aminoAcids = cat(3, NaN(20, 1), tmp(m.aminoAcidIndexs(1:20), 1, :));

            enzymeGene = cat(3, NaN(size(geneIdxs)), max(1, permute(MacromoleculeUtil.extractCopyNumberTimeCourse(c, states.Chromosome.polymerizedRegions, ...
                g.startCoordinates(geneIdxs), g.strands(geneIdxs), g.lengths(geneIdxs)), [1 3 2 4])));
            enzymeRNA = cat(3, NaN(size(matureRnaIdxs)), states.Rna.counts);
            enzymeMonomer = cat(3, NaN(size(monomerIdxs)), states.ProteinMonomer.counts(1:numel(monomerIdxs), 1, :, :));
            enzymeComplex = cat(3, NaN(size(complexIdxs)), states.ProteinComplex.counts(1:numel(complexIdxs), 1, :, :));
            reactionFlux = cat(3, NaN, states.MetabolicReaction.fluxs * 3600 * 1e-6);

            lk0 = full(sum(sum(states.Chromosome.polymerizedRegions, 2), 1)) / c.relaxedBasesPerTurn;
            superhelicalDensity = cat(3, NaN, (full(sum(sum(states.Chromosome.linkingNumbers, 2), 1)) - lk0) ./ lk0);
            %gyrase = ...
            %    + states.ProteinComplex.counts(end/2-1, 1, :, :) ...
            %    + states.ProteinComplex.counts(end-1, 1, :, :);
            %topoI = ...
            %    + states.ProteinMonomer.counts(end/2, 1, :, :) ...
            %    + states.ProteinMonomer.counts(end, 1, :, :);
            %topoIV = ...
            %    + states.ProteinComplex.counts(end/2, 1, :, :) ...
            %    + states.ProteinComplex.counts(end, 1, :, :);

            clear tmp states lk0;

            %% collect data
            figData = struct;
            figData.time = time;
            figData.mass = mass;
            figData.rnas = rnas;
            figData.proteins = proteins;
            figData.dntps = dntps;
            figData.ntps = ntps;
            figData.aminoAcids = aminoAcids;
            figData.growth = growth;
            figData.reactionFlux = reactionFlux;
            figData.enzymeComplex = enzymeComplex;
            figData.enzymeMonomer = enzymeMonomer;
            figData.enzymeRNA = enzymeRNA;
            figData.enzymeGene = enzymeGene;
            figData.ploidy = ploidy;
            figData.superhelicalDensity = superhelicalDensity;
        end

        function figData = mass(sim, simBatchDir, simNum, figHandle, normalized)
            import edu.stanford.covert.cell.sim.util.PlotUtil;
            import edu.stanford.covert.cell.sim.util.SimulationEnsemble;

            comp = sim.compartment;

            %% get data
            stateNames = {
                'Time' 'values'       ':' ':'
                'Mass' 'waterWt'      ':' comp.cellularIndexs(:)'
                'Mass' 'dnaWt'        ':' '-sum'
                'Mass' 'rnaWt'        ':' '-sum'
                'Mass' 'proteinWt'    ':' '-sum'
                'Mass' 'metaboliteWt' ':' [comp.cytosolIndexs; comp.membraneIndexs]
                };
            states = SimulationEnsemble.load(simBatchDir, stateNames, [], [], 1, 'extract', simNum);

            %% layout figure, plot data
            clf(figHandle);

            options = struct();
            options.titleStr = 'Cell Mass';

            options.plotFunc = @PlotUtil.plotStackedArea;
            axesHandle = PlotUtil.multiElementPlot(figHandle, 50, states.Time.values([1 end]), options);

            ydata = [
                permute(states.Mass.waterWt, [1 3 2])
                permute(states.Mass.proteinWt, [1 3 2])
                permute(states.Mass.dnaWt, [1 3 2])
                permute(states.Mass.metaboliteWt(:, 2, :), [1 3 2])
                permute(states.Mass.metaboliteWt(:, 1, :), [1 3 2])
                permute(states.Mass.rnaWt, [1 3 2])
                ] * 1e15;
            if nargin < 5 || ~normalized
                h = PlotUtil.plotStackedArea(axesHandle, permute(states.Time.values, [1 3 2]), ydata);
                ylabel(axesHandle, 'Mass (fg)');
                ylim(axesHandle, [0 max(sum(ydata, 1))]);
            else
                h = PlotUtil.plotNormalizedStackedArea(axesHandle, permute(states.Time.values, [1 3 2]), ydata);
                ylabel(axesHandle, 'Percent Mass');
                ylim(axesHandle, [0 100]);
            end
            legend(h, {'Water', 'Protein', 'DNA', 'Membrane', 'Metabolite', 'RNA'}, 'Location', 'SouthEast');

            PlotUtil.offsetYAxes(axesHandle, 0.02);

            clear ydata;

            %% collect data
            figData = struct;
            figData.time = states.Time.values;
            figData.waterWt = states.Mass.waterWt;
            figData.proteinWt = states.Mass.proteinWt;
            figData.dnaWt = states.Mass.dnaWt;
            figData.metaboliteWt = states.Mass.metaboliteWt;
            figData.rnaWt = states.Mass.rnaWt;
        end

        function figData = massNormalized(sim, simBatchDir, simNum, figHandle)
            import edu.stanford.covert.cell.sim.analysis.SingleCell;

            figData = SingleCell.mass(sim, simBatchDir, simNum, figHandle, true);
        end

        function figData = metabolites(sim, simBatchDir, simNum, figHandle, normalized)
            import edu.stanford.covert.cell.sim.util.PlotUtil;
            import edu.stanford.covert.cell.sim.util.SimulationEnsemble;
            import edu.stanford.covert.util.ConstantUtil;

            comp = sim.compartment;
            m = sim.state('Metabolite');
            metIdxs = [
                m.ntpIndexs
                m.ndpIndexs
                m.nmpIndexs
                m.dntpIndexs
                m.aminoAcidIndexs
                m.phosphateIndexs
                m.diphosphateIndexs
                m.hydrogenIndexs
                ];

            %% get data
            stateNames = {
                'Time'       'values'  ':'      ':'
                'Geometry'   'volume'  ':'      ':'
                'Metabolite' 'counts'  metIdxs  comp.cytosolIndexs
                };
            states = SimulationEnsemble.load(simBatchDir, stateNames, [], [], 1, 'extract', simNum);
            time = permute(states.Time.values, [1 3 2]) / 3600;

            states.Metabolite.counts = ...
                states.Metabolite.counts ./ ...
                states.Geometry.volume(ones(size(metIdxs)), :, :) / ...
                ConstantUtil.nAvogadro * 1000;

            ntps = permute(states.Metabolite.counts(1:4, :, :), [1 3 2]);
            ndps = permute(states.Metabolite.counts(5:8, :, :), [1 3 2]);
            nmps = permute(states.Metabolite.counts(9:12, :, :), [1 3 2]);
            dntps = permute(states.Metabolite.counts(13:16, :, :), [1 3 2]);
            aas = permute(states.Metabolite.counts(17:37, :, :), [1 3 2]);
            others = permute(states.Metabolite.counts(38:40, :, :), [1 3 2]);

            %% layout figure
            clf(figHandle);

            axesHandles = PlotUtil.multiElementPlot(figHandle, {10 * ones(3, 1) 10 * ones(3, 1)}, time([1 end]));

            %% plot data
            if nargin >= 5 && normalized
                plotFcn = @PlotUtil.plotNormalizedStackedArea;
            else
                plotFcn = @PlotUtil.plotStackedArea;
            end

            %ntp
            axesHandle = axesHandles{1}(1);
            h = plotFcn(axesHandle, time, ntps);
            legend(h, m.wholeCellModelIDs(m.ntpIndexs), 'Location', 'SouthEast');

            %ndp
            axesHandle = axesHandles{1}(2);
            h = plotFcn(axesHandle, time, ndps);
            legend(h, m.wholeCellModelIDs(m.ndpIndexs), 'Location', 'SouthEast');

            %nmp
            axesHandle = axesHandles{1}(3);
            h = plotFcn(axesHandle, time, nmps);
            legend(h, m.wholeCellModelIDs(m.nmpIndexs), 'Location', 'SouthEast');

            %dNTP
            axesHandle = axesHandles{2}(1);
            h = plotFcn(axesHandle, time, dntps);
            legend(h, m.wholeCellModelIDs(m.dntpIndexs), 'Location', 'SouthEast');

            %AA
            axesHandle = axesHandles{2}(2);
            plotFcn(axesHandle, time, aas);

            %other
            axesHandle = axesHandles{2}(3);
            h = plotFcn(axesHandle, time, others);
            legend(h, m.wholeCellModelIDs([m.phosphateIndexs; m.diphosphateIndexs; m.hydrogenIndexs]), 'Location', 'SouthEast');

            %labels
            for i = 1:numel(axesHandles)
                for j = 1:numel(axesHandles{i})
                    if nargin >= 5 && normalized
                        ylabel(axesHandles{i}(j), {'Relative' 'Copy Number'});
                        ylim(axesHandles{i}(j), [0 100]);
                    else
                        ylabel(axesHandles{i}(j), 'Concentration (mM)');
                    end
                end
            end

            %% ylabels
            PlotUtil.alignYAxesLabels(axesHandles);
            PlotUtil.offsetYAxes(axesHandles, 0.02);

            %% collect data
            figData = struct;
            figData.time = time;
            figData.ntps = ntps;
            figData.ndps = ndps;
            figData.nmps = nmps;
            figData.dntps = dntps;
            figData.aas = aas;
            figData.others = others;
        end

        function figData = metabolitesNormalized(sim, simBatchDir, simNum, figHandle)
            import edu.stanford.covert.cell.sim.analysis.SingleCell;

            figData = SingleCell.metabolites(sim, simBatchDir, simNum, figHandle, true);
        end

        function figData = ribosomes(sim, simBatchDir, simNum, figHandle, useGenomicCoordinates)
            import edu.stanford.covert.cell.sim.util.PlotUtil;
            import edu.stanford.covert.cell.sim.util.SimulationEnsemble;

            %get constants
            g = sim.gene;
            c = sim.state('Chromosome');
            rib = sim.state('Ribosome');

            if nargin >= 5 && useGenomicCoordinates
                startCoordinates = g.startCoordinates(g.mRNAIndexs);
                lengths = g.lengths(g.mRNAIndexs);
                strands = g.strands(g.mRNAIndexs);
            else
                startCoordinates = cumsum(g.lengths(g.mRNAIndexs)+50);
                lengths = g.lengths(g.mRNAIndexs);
                strands = ones(size(g.mRNAIndexs));
            end

            %get data
            stateNames = {
                'Time'      'values'
                'Ribosome'  'states'
                'Ribosome'  'boundMRNAs'
                'Ribosome'  'mRNAPositions'
                };
            states = SimulationEnsemble.load(simBatchDir, stateNames, [], [], 1, 'extract', simNum);
            time = permute(states.Time.values, [1 3 2]) / 3600;
            ribStates = permute(states.Ribosome.states, [1 3 2]);
            boundMRNAs = permute(states.Ribosome.boundMRNAs, [1 3 2]);
            mRNAPos = permute(states.Ribosome.mRNAPositions, [1 3 2]);

            %layout figure
            clf(figHandle);
            options = struct;
            axesHandle = PlotUtil.multiElementPlot(figHandle, 30, time([1 end]), options);

            %plot data
            colorOrder = PlotUtil.getRedGreenColorOrder(randperm(size(ribStates, 1)));
            for i = 1:size(ribStates, 1)
                pos = NaN(size(time));
                tfs = ribStates(i, :) == rib.activeValue & mRNAPos(i, :) > 0;
                pos(tfs) = ...
                    + startCoordinates(boundMRNAs(i, tfs)) ...
                    + (strands(boundMRNAs(i, tfs)) == 2) .* (lengths(boundMRNAs(i, tfs)) - 1) ...
                    + (2 * (2 - strands(boundMRNAs(i, tfs))) - 1) .* (3 * mRNAPos(i, tfs)' - 1);
                h = plot(axesHandle, time, pos);
                set(h, 'Color', colorOrder(i, :));
            end

            if nargin >= 5 && useGenomicCoordinates
                ylabel(axesHandle, 'Ribosome Position');
                ylim(axesHandle, [1 c.sequenceLen]);
                set(axesHandle, ...
                    'YTick', [1 c.terCPosition c.sequenceLen], ...
                    'YTickLabel', {'OriC' 'TerC' 'Oric'});
                PlotUtil.offsetYAxes(axesHandle, 0.02);
            else
                tmp = xlim(axesHandle);
                x = tmp(1) - 0.02 * range(tmp);
                xlim(axesHandle, [x tmp(2)]);
                ylim(axesHandle, [min(startCoordinates) max(startCoordinates + lengths - 1)]);
                line(x * ones(3 * numel(startCoordinates), 1), ...
                    reshape([startCoordinates  startCoordinates + lengths - 1  NaN(size(startCoordinates))]', [], 1), ...
                    'Parent', axesHandle, 'Color', 'k')
                set(axesHandle, 'YTick', [], 'visible', 'off');
                text(tmp(1) - 0.03 * range(tmp), mean(ylim(axesHandle)), 'Ribosome Position', ...
                    'Parent', axesHandle, ...
                    'rotation', 90, ...
                    'HorizontalAlignment', 'center', ...
                    'VerticalAlignment', 'bottom');
            end

            %% collect data
            figData = struct;
            figData.time = time;
            figData.ribStates = ribStates;
            figData.boundMRNAs = boundMRNAs;
            figData.mRNAPos = mRNAPos;
        end

        function figData = ribosomesGenomicCoordinates(sim, simBatchDir, simNum, figHandle)
            import edu.stanford.covert.cell.sim.analysis.SingleCell;

            figData = SingleCell.ribosomes(sim, simBatchDir, simNum, figHandle, true);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
%% Image Moments
% We learn to calculate the moments of an image.
%
% In this sample you will learn how to use the following OpenCV functions:
%
% * <matlab:doc('cv.moments') cv.moments>
% * <matlab:doc('cv.contourArea') cv.contourArea>
% * <matlab:doc('cv.arcLength') cv.arcLength>
%
% Sources:
%
% * <https://docs.opencv.org/3.4.0/d0/d49/tutorial_moments.html>
% * <https://github.com/opencv/opencv/blob/3.4.0/samples/cpp/tutorial_code/ShapeDescriptors/moments_demo.cpp>
%

function varargout = moments_demo_gui(im)
    % load source image
    if nargin < 1
        src = imread(fullfile(mexopencv.root(),'test','monster.jpg'));
    elseif ischar(im)
        src = imread(im);
    else
        src = im;
    end

    % Convert image to gray and blur it
    if size(src,3) == 3
        src = cv.cvtColor(src, 'RGB2GRAY');
    end
    src = cv.blur(src, 'KSize',[3 3]);

    % create the UI
    h = buildGUI(src);
    if nargout > 0, varargout{1} = h; end
end

function onChange(~,~,h)
    %ONCHANGE  Event handler for UI controls

    % retrieve current values from UI controls
    thresh = round(get(h.slid, 'Value'));
    set(h.txt, 'String',sprintf('Canny thresh: %3d',thresh));

    % Detect edges using Canny
    canny_output = cv.Canny(h.src, [thresh thresh*2], 'ApertureSize',3);

    % Find contours
    contours = cv.findContours(canny_output, 'Mode','Tree', 'Method','Simple');

    % Get the moments and compute the mass center
    mu = cell(size(contours));
    mc = cell(size(contours));
    for i=1:numel(contours)
        mu{i} = cv.moments(contours{i}, 'BinaryImage',false);
        mc{i} = [mu{i}.m10, mu{i}.m01] ./ mu{i}.m00;
    end

    % Draw contours
    drawing = zeros([size(canny_output) 3], 'uint8');
    for i=1:numel(contours)
        clr = randi([0 255], [1 3], 'uint8');
        drawing = cv.drawContours(drawing, contours, 'ContourIdx',i-1, ...
            'Color',clr, 'Thickness',2);
        drawing = cv.circle(drawing, mc{i}, 4, ...
            'Color',clr, 'Thickness','Filled');
    end

    % show result
    set(h.img, 'CData',drawing);
    drawnow;

    % Calculate the area with the moments 00 and compare with the result of
    % the OpenCV function
    fprintf('\t Info: Area and Contour Length \n');
    for i=1:numel(contours)
        fprintf([' * Contour[%d] - Area (M_00) = %.2f - Area OpenCV: %.2f' ...
            ' - Length: %.2f \n'], i, mu{i}.m00, ...
            cv.contourArea(contours{i}), ...
            cv.arcLength(contours{i}, 'Closed',true));
    end
end

function h = buildGUI(img)
    %BUILDGUI  Creates the UI

    % parameters
    thresh = 100;
    max_thresh = 255;
    sz = size(img);
    sz(2) = max(sz(2), 250);  % minimum figure width

    % build the user interface (no resizing to keep it simple)
    h = struct();
    h.src = img;
    h.fig = figure('Name','Image Moments Demo', ...
        'NumberTitle','off', 'Menubar','none', 'Resize','off', ...
        'Position',[200 200 sz(2) sz(1)+29]);
    if ~mexopencv.isOctave()
        %HACK: not implemented in Octave
        movegui(h.fig, 'center');
    end
    h.ax = axes('Parent',h.fig, 'Units','pixels', 'Position',[1 30 sz(2) sz(1)]);
    if ~mexopencv.isOctave()
        h.img = imshow(img, 'Parent',h.ax);
    else
        %HACK: https://savannah.gnu.org/bugs/index.php?45473
        axes(h.ax);
        h.img = imshow(img);
    end
    h.txt = uicontrol('Parent',h.fig, 'Style','text', 'FontSize',11, ...
        'Position',[5 5 130 20], 'String',sprintf('Canny thresh: %3d',thresh));
    h.slid = uicontrol('Parent',h.fig, 'Style','slider', 'Value',thresh, ...
        'Min',0, 'Max',max_thresh, 'SliderStep',[1 10]./(max_thresh-0), ...
        'Position',[135 5 sz(2)-135-5 20]);

    % hook event handlers, and trigger default start
    set(h.slid, 'Callback',{@onChange,h}, ...
        'Interruptible','off', 'BusyAction','cancel');
    onChange([],[],h);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [ CostMatric,Cost ] =  MinimalDirectedMSF( CostMatric )
% MST on a directed graph
% Chu-Liu/Edmonds Algorithm:
%1 Discard the arcs entering the root if any; For each node other than the root, select the entering arc with the highest cost; Let the selected n-1 arcs be the set S.
%  If no cycle formed, G( N,S ) is a MST. Otherwise, continue.
%2 For each cycle formed, contract the nodes in the cycle into a pseudo-node (k), and modify the cost of each arc which enters a node (j) in the cycle from some node (i)
%  outside the cycle according to the following equation.  c(i,k)=c(i,j)-(c(x(j),j)-min_{j}(c(x(j),j)) where c(x(j),j) is the cost of the arc in the cycle which enters j.
%3 For each pseudo-node, select the entering arc which has the smallest modified cost; Replace the arc which enters the same real node in S by the new selected arc.
%4 Go to step 2 with the contracted graph.

% Pay attention, we extend the Chu-Liu algorithm into considering
% forest.And we consider of discarding the premise of giving the root here.
% CostMatric must be symmetric

  [ Number, Component ] = conncomp(  biograph( CostMatric ),'Weak', true );
  for k = 1:Number
      LocalNodes = find( Component == k );
      Dim = size( LocalNodes,2 );
      if Dim == 1
         continue
      elseif Dim == 2
          if CostMatric( LocalNodes(1),LocalNodes(2) ) > CostMatric( LocalNodes(2),LocalNodes(1) )
              CostMatric( LocalNodes(2),LocalNodes(1) ) = 0;
          else
              CostMatric( LocalNodes(1),LocalNodes(2) ) = 0;
          end
      elseif Dim > 2
      ComponentCost = CostMatric( LocalNodes,LocalNodes );
      MinTree = zeros( Dim ); MinTreeCost = Inf;
      for m = 1:Dim
          Root = m;
          LocalCostMatric = ComponentCost;
          LocalCostMatric( :,Root ) = zeros( Dim,1 );
          [ LocalTree,LocalCostTree ] =  DirectedMinimalSpanningTree( LocalCostMatric,Root );
          if MinTreeCost > LocalCostTree
             MinTree = LocalTree;
             MinTreeCost = LocalCostTree ;
          end
      end
      CostMatric( LocalNodes,LocalNodes ) = MinTree;
      end
   end
   Cost = sum( sum( CostMatric ));
end


text:
% Gives the spectrum of a derivative effect (i.e. a zero at zero frequency)
%
% Input
%  dftlen : the number of bin in spectrum (full DFT length).
%  fs     : [Hz] The sampling frequency
%
% Output
%  S      : The spectrum of the derivative effect
%
% Copyright (c) 2011 University of Crete - Computer Science Department
%
% License
%  This file is under the LGPL license,  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. This file 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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function S = spec_derivative(dftlen, fs)

    % The half-spectrum of a zero at frequency zero is:
    hS = fs*(2i*pi/dftlen)*(0:dftlen/2).';

    if mod(dftlen,2)==1
        % If DFT length is odd
        S = [hS; hS(end:-1:2)];
    else
        % If DFT length is even
        hS(end) = 0;
        S = hspec2spec(hS);
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [patches, pdm, clmParams, early_term_params] = Load_CECLM_general()
%LOAD_PATCH_EXPERTS Summary of this function goes here
%   Detailed explanation goes here

    % Load the patch experts/local detectors
    [patches] = Load_CECLM_Patch_Experts( '../models/cen/', 'cen_patches_*_general.mat');

    % the default PDM to use
    pdmLoc = ['../models/pdm/pdm_68_aligned_wild.mat'];

    load(pdmLoc);

    pdm = struct;
    pdm.M = double(M);
    pdm.E = double(E);
    pdm.V = double(V);

    clmParams = struct;

    clmParams.window_size = [25,25; 23,23; 21,21; 21,21];
    clmParams.numPatchIters = size(clmParams.window_size,1);

    clmParams.regFactor = 0.9*[35, 27, 20, 20];
    clmParams.sigmaMeanShift = 1.5*[1.25, 1.375, 1.5, 1.5];
    clmParams.tikhonov_factor = [2.5, 5, 7.5, 7.5];

    clmParams.startScale = 1;
    clmParams.num_RLMS_iter = 10;
    clmParams.fTol = 0.01;
    clmParams.useMultiScale = true;
    clmParams.use_multi_modal = 1;
    clmParams.multi_modal_types  = patches(1).multi_modal_types;
    clmParams.numPatchIters = 4;

    % As the orientations are not equally reliable reweigh them
    load('../models/cen/cen_general_mapping.mat');
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function [X, A] = PCA_stochastic(A, k)
% Example of stochastic gradient algorithm in Manopt on a PCA problem.
%
% PCA (principal component analysis) on a dataset A of size nxd consists
% in solving
%
%   minimize_X  f(X) = -.5*norm(A*X, 'fro')^2 / n,
%
% where X is a matrix of dimension dxk with orthonormal columns. This
% is equivalent to finding k dominant singular vectors of A, or k top
% eigenvectors of A'*A.
%
% If n is large, this computation can be expensive. Thus,  stochastic
% gradient algorithms take the point of view that f(X) is a sum of many (n)
% terms: each term involves only one of the n rows of A.
%
% To make progress, it may be sufficient to optimize with respect to a
% subset of the terms at each iteration. This way, each individual
% iteration can be very cheap. In particular, individual operations have
% cost independent of n, because f or its gradient need never be evaluated
% completely (or at all in the case of f.)
%
% Stochastic gradient algorithms (this implementation in particular) are
% sensitive to proper parameter tuning. See in code.

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Bamdev Mishra and Nicolas Boumal, Sept. 6, 2017
% Contributors:
%
% Change log:
%

    % If none is given, generate a random data set: n samples in R^d
    if ~exist('A', 'var') || isempty(A)
        d = 1000;
        n = 100000;
        fprintf('Generating data...');
        A = randn(n, d)*diag([[15 10 5], ones(1, d-3)]);
        fprintf(' done (size: %d x %d).\n', size(A));
    else
        [n, d] = size(A);
    end

    % Pick a number of component to compute
    if ~exist('k', 'var') || isempty(k)
        k = 3;
    end

    % We are looking for k orthonormal vectors in R^d: Stiefel manifold.
    problem.M = stiefelfactory(d, k);

    % The cost function to minimize is a sum of n terms. This parameter
    % must be set for stochastic algorithms.
    problem.ncostterms = n;

    % We do not need to specify how to compute the value of the cost
    % function (stochastic algorithms never use this). All we need is to
    % specify how to compute the gradient of the cost function, where the
    % sum is restricted to a subset of the terms (a sample). Notice that we
    % specify a partial Euclidean gradient (hence the 'e' in partialegrad).
    % This way, Manopt will automatically convert the Euclidean vector into
    % a proper Riemannian partial gradient, in the tangent space at X.
    % In particular, if sample = 1:n, then the partial gradient corresponds
    % to the actual (complete) gradient.
    problem.partialegrad = @partialegrad;
    function G = partialegrad(X, sample)

        % X is an orthonormal matrix of size dxk
        % sample is a vector if indices between 1 and n: a subset
        % Extract a subset of the dataset
        Asample = A(sample, :);

        % Compute the gradient of f restricted to that sample
        G = -Asample'*(Asample*X);
        G = G / n;

    end

    % If one wants to use checkgradient to verify one's work, then it is
    % necessary to specify the cost function as well, as below.
    % problem.cost = @(X) -.5*norm(A*X, 'fro')^2 / n;
    % checkgradient(problem); pause;

    % To have the solver record statistics every x iterations, set
    % options.checkperiod to x. This will record simple quantities which
    % are almost free to compute (namely, elapsed time and step size of the
    % last step.) To record more sophisticated quantities, you can use
    % options.statsfun as usual. Time spent computing these statistics is
    % not counted in times reported in the info structure returned by the
    % solver.
    options.checkperiod = 10;
    options.statsfun = statsfunhelper('metric', @(X) norm(A*X, 'fro'));

    % Set the parameters for the solver: stochastic gradient algorithms
    % tend to be quite sensitive to proper tuning, especially regarding
    % step size selection. See the solver's documentation for details.
    options.maxiter = 200;
    options.batchsize = 10;
    % options.stepsize_type = 'decay';
    options.stepsize_init = 1e2;
    options.stepsize_lambda = 1e-3;
    options.verbosity = 2;

    % Run the solver
    [X, info] = stochasticgradient(problem, [], options);

    % Plot the special metric recorded by options.statsfun
    plot([info.iter], [info.metric], '.-');
    xlabel('Iteration #');
    ylabel('Frobenius norm of A*X');
    title('Convergence of stochasticgradient on stiefelfactory for PCA');

    % Add to that plot a reference: the globally optimal value attained if
    % the true dominant singular vectors are computed.
    fprintf('Running svds... ');
    t = tic();
    [V, ~] = svds(A', k);
    fprintf('done: %g [s] (note: svd may be faster)\n', toc(t));
    hold all;
    bound = norm(A*V, 'fro');
    plot([info.iter], bound*ones(size([info.iter])), '--');
    hold off;

    legend('Algorithm', 'SVD bound', 'Location', 'SouthEast');

end

author: NicolasBoumal

repo: manopt