text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = INVERSEKINEMATIC_FLEXION_N2(robot, T)
%   Solves the inverse kinematic problem for the ABB IRB 140 robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=INVERSEKINEMATIC_FLEXION_N2 returns 8 possible solutions, thus,
%   Q is a 6x8 matrix where each column stores 6 feasible joint values.
%
%
%   Example code:
%
%   >>epson=load_robot('EPSON', 'FLEXION_N2');
%   >>q = [0 0 0 0 0 0];
%   >>T = directkinematic(epson, q);
%
%   %Call the inversekinematic for this robot
%
%   >> qinv = inversekinematic(epson, T);
%
%   check that all of them are feasible solutions!
%   and every Ti equals T
%
%   for i=1:8,
%        Ti = directkinematic(epson, qinv(:,i))
%   end
%
%	See also DIRECTKINEMATIC.
%
%   Author: Arturo Gil Aparicio
%           Universitas Miguel Hernandez, SPAIN.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE 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.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function q = inversekinematic_flexion_n2(robot, T)

Ts=[-1 0 0 0; 0 -1 0 0; 0 0 -1 0.9; 0 0 0 1];
T=inv(Ts)*T;
%initialize q,
%eight possible solutions are generally feasible
q=zeros(6,8);

%Evaluate the parameters
d = eval(robot.DH.d);

%See geometry at the reference for this robot
L6=d(6);

%A1 = a(1);

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Compute the position of the wrist, being W the Z component of the end effector's system
W = T(1:3,3);

% Pm: wrist position
Pm = [Px Py Pz]' - L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution
q1=atan2(Pm(2), Pm(1))+pi/2;

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);

q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);

q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%Arrange solutions, there are 8 possible solutions so far.
% if q1 is a solution, q1* = q1 + pi is also a solution.
% For each (q1, q1*) there are two possible solutions
% for q2 and q3 (namely, elbow up and elbow up solutions)
% So far, we have 4 possible solutions. Howefer, for each triplet (theta1, theta2, theta3),
% there exist two more possible solutions for the last three joints, generally
% called wrist up and wrist down solutions. For this reason,
%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];

%leave only the real part of the solutions
q=real(q);

%Note that in this robot, the joint q3 has a non-simmetrical range. In this
%case, the joint ranges from 60 deg to -219 deg, thus, the typical normalizing
%step is avoided in this angle (the next line is commented). When solving
%for the orientation, the solutions are normalized to the [-pi, pi] range
%only for the theta4, theta5 and theta6 joints.

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));

% solve for the last three joints
% for any of the possible combinations (theta1, theta2, theta3)
for i=1:2:size(q,2),
    % use solve_spherical_wrist2 for the particular orientation
    % of the systems in this EPSON robot
    % use either the geometric or algebraic method.
    % the function solve_spherical_wrist2 is used due to the relative
    % orientation of the last three DH reference systems.

    %use either one algebraic method or the geometric
    qtemp = solve_spherical_wrist_flexion_n2(robot, q(:,i), T, 1, 'geometric'); %wrist up
    %qtemp = solve_spherical_wrist_flexion_n2(robot, q(:,i), T, 1,'algebraic'); %wrist up
    %qtemp = solve_spherical_wrist(robot, q(:,i), T, 1, 'geometric'); %wrist up
    %qtemp = solve_spherical_wrist(robot, q(:,i), T, 1,'algebraic'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    qtemp = solve_spherical_wrist_flexion_n2(robot, q(:,i), T, -1, 'geometric'); %wrist down
    %qtemp = solve_spherical_wrist_flexion_n2(robot, q(:,i), T, -1,'algebraic'); %wrist up
    %qtemp = solve_spherical_wrist(robot, q(:,i), T, -1, 'algebraic'); %wrist down
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
d = eval(robot.DH.d);
a = eval(robot.DH.a);

%See geometry
L2=a(2);
L3=d(4);

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = atan2(-p1(2), p1(1));
gamma = (acos((L2^2+r^2-L3^2)/(2*r*L2)));

if ~isreal(gamma)
    disp('WARNING:inversekinematic_flexion_n2: the point is not reachable for this configuration, imaginary solutions');
    %gamma = real(gamma);
end

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta3
q2(1) = pi/2 - beta - gamma; %elbow up
q2(2) = pi/2 - beta + gamma; %elbow down

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
d = eval(robot.DH.d);
a = eval(robot.DH.a);

%See geometry
L2=a(2);
L3=d(4);

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

eta = (acos((L2^2 + L3^2 - r^2)/(2*L2*L3)));

if ~isreal(eta)
   disp('WARNING:inversekinematic_flexion_n2: the point is not reachable for this configuration, imaginary solutions');
   %eta = real(eta);
end

%return two possible solutions
%elbow up and elbow down solutions
%the order here is important
q3(1) = pi - eta;
q3(2) = eta - 3*pi/2 + pi/2;

repo: ARTE


text:
% BARH Horizontal bar graph.
%    BARH(X,Y) draws the columns of the M-by-N matrix Y as M groups of
%    N horizontal bars. The vector X must not have duplicate values.
%
%    BARH(Y) uses the default value of X=1:M.  For vector inputs,
%    BARH(X,Y) or BARH(Y) draws LENGTH(Y) bars.  The colors are set by
%    the colormap.
%
%    BARH(X,Y,WIDTH) or BARH(Y,WIDTH) specifies the width of the
%    bars. Values of WIDTH > 1, produce overlapped bars.  The
%    default value is WIDTH=0.8.
%
%    BARH(...,'grouped') produces the default horizontal grouped bar chart.
%    BARH(...,'stacked') produces a horizontal stacked bar chart.
%    BARH(...,LINESPEC) uses the line color specified (one of 'rgbymckw').
%
%    BARH(AX,...) plots into AX instead of GCA.
%
%    H = BARH(...) returns a vector of handles to barseries objects.
%
%    Use SHADING FACETED to put edges on the bars.  Use SHADING FLAT to
%    turn them off.
%
%    Examples: subplot(3,1,1), barh(rand(10,5),'stacked'), colormap(cool)
%              subplot(3,1,2), barh(0:.25:1,rand(5),1)
%              subplot(3,1,3), barh(rand(2,3),.75,'grouped')
%
%    See also PLOT, BAR, BAR3, BAR3H, HISTOGRAM.
%
%    Copyright 1984-2017 The MathWorks, Inc.
%
%    Reference page in Doc Center
%       doc barh
%
%    Other functions named barh
%
%       fints/barh    ts/barh
%

author: jmaih


text:
function [bval] = f2b(fval,bounds,bits)
% function [bval] = f2b(fval,bounds,bits)
%
% Return the binary representation of the float number fval.
%
% fval   - the float representation of the number
% bval   - the binary representation of the number
% bounds - the bounds on the variables
% bits   - the number of bits to represent each variable

% Binary and Real-Valued Simulation Evolution for Matlab
% Copyright (C) 1996 C.R. Houck, J.A. Joines, M.G. Kay
%
% C.R. Houck, J.Joines, and M.Kay. A genetic algorithm for function
% optimization: A Matlab implementation. ACM Transactions on Mathmatical
% Software, Submitted 1996.
%
% 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 1, 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. A copy of the GNU
% General Public License can be obtained from the
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

scale=(2.^bits-1)./ (bounds(:,2)-bounds(:,1))'; %The range of the variables
numV=size(bounds,1);
cs=[0 cumsum(bits)];
bval=[];
for i=1:numV
  fval(i)=(fval(i)-bounds(i,1)) * scale(i);
  bval=[bval rem(floor(fval(i)*pow2(1-bits(i):0)),2)];
end

author: HuangCongQing


text:
% Script demo_copy_data
% Shows functionality of MrCopyData for deep cloning and recursive operations
%
%  demo_copy_data
%
% Since MrCopyData has no properties, we created a dummy class
% MrCopyDataTest with the following properties:
%   - numeric/string values, cells of nums/strings,
%   - MrCopyDataTest values and
%   - cells of MrCopyDataTest
% to allow for a minimum comparison case of the recursive operations
% provided by MrCopyData. The MrCopyDataTest entries are recursive objects,
% but with one recursion depth less, i.e. at some point in the hierarchy,
% copydata_value and _cell will become empty
%
%   See also MrCopyData MrCopyDataTest

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2017-07-20
% Copyright (C) 2017 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

clear;
close all;
clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Create a deep MrCopyDataObject of given recursion depth
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nRecursionDepth = 2;
Y = MrCopyDataTest(nRecursionDepth);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Perform some unary operations for fun
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Y.print
Y.find('MrCopyDataTest'); % recursive find, all props of this class
Y.get_nonempty_fields;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Clone object and check comparison/diff works
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Z = Y.copyobj();
isequal(Y,Z) % should be the same...
dYZ = Y.diffobj(Z); % should have all-empty values...but resembles the recursive prop structure of the original object

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Manipulate clone and check comparison/diff works
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
Z.copydata_value.clear(); % now sth is missing...
[dYZ, isYZEqual] = Y.diffobj(Z); %  "Y - Z", per property, chooses Y's prop where Z differs

Y.print_diff(Z);

isequal(Y,Z) % should be false...
isequal(Z,Y) % should be false...

% update all non-empty properties from other object
Z.update_properties_from(dYZ,1);
isequal(Y,Z) % TODO: bugfix! should be true again...

author: translationalneuromodeling

repo: tapas


text:
function [gradp,varlist] = gradient(sp)
% sympoly/gradient: gradient vector of a sympoly
% usage: [gradp,varlist] = gradient(sp);
%
% arguments: (input)
%     sp - scalar sympoly object
%
% arguments: (output)
%  gradp - sympoly (row vector) object containing the gradient vector
%
%  varlist - cell array of variable names used for the gradient

if numel(sp)>1
  error 'Gradient only works for scalar sympoly objects.'
end

varlist = setdiff(sp.Var,{''});
nvar=length(varlist);

if nvar==0
  % sympoly had no variables, i.e., it was a constant
  gradp=sympoly(0);
else
  % loop over the variables
  gradp = sympoly(zeros(1,nvar));
  for i=1:nvar
    gradp(i)=diff(sp,1,varlist{i});
  end
end

author: Sable

repo: mcbench-benchmarks


text:
% mtimesx_sparse does sparse matrix multiply of two inputs
%******************************************************************************
%
%  MATLAB (R) is a trademark of The Mathworks (R) Corporation
%
%  Function:    mtimesx_sparse
%  Filename:    mtimesx_sparse.m
%  Programmer:  James Tursa
%  Version:     1.00
%  Date:        September 27, 2009
%  Copyright:   (c) 2009 by James Tursa, All Rights Reserved
%
%  This code uses the BSD License:
%
%  Redistribution and use in source and binary forms, with or without
%  modification, are permitted provided that the following conditions are
%  met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
%  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%  POSSIBILITY OF SUCH DAMAGE.
%
%--
%
% mtimesx_sparse is a helper function for mtimesx and is not intended to be called
% directly by the user.
%
% ---------------------------------------------------------------------------------------------------------------------------------

function result = mtimesx_sparse(a,transa,b,transb)
if( transa == 'N' )
    if( transb == 'N' )
        result = a * b;
    elseif( transb == 'G' )
        result = a * conj(b);
    elseif( transb == 'T' )
        result = a * b.';
    else
        result = a * b';
    end
elseif( transa == 'G' )
    if( transb == 'N' )
        result = conj(a) * b;
    elseif( transb == 'G' )
        result = conj(a) * conj(b);
    elseif( transb == 'T' )
        result = conj(a) * b.';
    else
        result = conj(a) * b';
    end
elseif( transa == 'T' )
    if( transb == 'N' )
        result = a.' * b;
    elseif( transb == 'G' )
        result = a.' * conj(b);
    elseif( transb == 'T' )
        result = a.' * b.';
    else
        result = a.' * b';
    end
else
    if( transb == 'N' )
        result = a' * b;
    elseif( transb == 'G' )
        result = a' * conj(b);
    elseif( transb == 'T' )
        result = a' * b.';
    else
        result = a' * b';
    end
end
end

author: martin-danelljan

repo: ECO


text:
function gdoc(varargin)

% function gdoc(functionName)
%------------------------------------------------------------------------
% View XML files using a brower embedded in a figure window or within
% MATLAB brower.
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2016/05/13 %Updated, renamed from febView to xmlView
% 2017/08/18 %Created varargin with viewer option. Create temp file to
% allow viewing of domNode.
%------------------------------------------------------------------------
%%

switch nargin
    case 1
        topic=varargin{1};
    case 0
        topic='';
end

if isempty(topic)
   web('GIBBON_product_page.html', '-helpbrowser'); %Display main GIBBON page
else
    documentationPath=fullfile(fileparts(fileparts(mfilename('fullpath'))),'docs','html');

    helpCheckName=fullfile(documentationPath,['HELP_',topic,'.html']);
    demoCheckName=fullfile(documentationPath,['DEMO_',topic,'.html']);

    if exist(helpCheckName,'file')
        web(helpCheckName, '-helpbrowser');
    elseif exist(demoCheckName,'file')
        web(demoCheckName, '-helpbrowser');
    else
        disp('No matching HELP or DEMO files found. Searching all documentation instead');
        docsearch(topic);
    end
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function restorefig(h,old)
    %RESTOREFIG  Restore a figure's properties
    %   RESTOREFIG(H,OLD) restores the properties of H specified in
    %   OLD. The state-difference structure OLD is the output of the
    %   function APPLYTOFIG. Deleting object handles after the call to
    %   APPLYFIG might cause undesired behavior.
    %
    %   See also EXPORTFIG, PREVIEWFIG, APPLYTOFIG.

    %  Copyright 2000 Ben Hinkle
    %  Email bug reports and comments to bhinkle@mathworks.com

    for n=1:length(old.objs)
        try
            if ~iscell(old.values{n})  &&  iscell(old.prop{n})
                old.values{n} = {old.values{n}};
            end
            set(old.objs{n}, old.prop{n}, old.values{n});
        end
    end

author: CelsoReyes


text:
function res = bnt_to_mat(data,misv)
% res = bnt_to_mat(data,misv)
%

if nargin <1, error('Requires at least 1 argument.'), end
if nargin == 1,	misv=-9999; end
taille=size(data);
long=taille(1);
larg=taille(2);
for i=1:long
  for j=1:larg
    if ~isempty(data{i,j})
      res(i,j)=data{i,j};
    else
      res(i,j)=misv;
    end
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
classdef UnfittedMesh_FlatSurface < UnfittedMesh_AbstractBuilder
    properties (GetAccess = public, SetAccess = private)
        unfittedType = 'INTERIOR';
        maxSubcells = 6;
        nnodesSubcell = 3;

        subcellsMesher      = SubcellsMesher_Interior;
        cutPointsCalculator = CutPointsCalculator;
        meshPlotter         = MeshPlotter_Interior_2D;
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%%*****************************************************************
%% schurmat_lblk: compute A*D*A'
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function [schur,UU,EE] = schurmat_lblk(blk,At,par,schur,UU,EE,p,dd)

   global idxdenAl

   iter = par.iter;
   n = sum(blk{p,2});

   if (iter==1)
      idxdenAl{p} = checkdense(At{p}');
   end
   ddsch = dd{p};
   if ~isempty(idxdenAl{p});
      idxden = idxdenAl{p};
      len = length(idxden);
      Ad = At{p}(idxden,:)' *spdiags(sqrt(ddsch(idxden)),0,len,len);
      UU = [UU, Ad];
      if isempty(EE)
         count = 0;
      else
         count = max(max(EE(:,1)),max(EE(:,2)));
      end
      tmp = count + [1:len]';
      EE = [EE; [tmp, tmp, -ones(len,1)] ];
      ddsch(idxden) = zeros(len,1);
   end
   schurtmp = At{p}' *spdiags(ddsch,0,n,n) *At{p};
   schur = schur + schurtmp;
%%*******************************************************************

author: shamilmamedov


text:
function h=hmatrix(p,xx,yy,dd,hh,varargin)

%   Copyright (C) 2004-2006 Per-Olof Persson. See COPYRIGHT.TXT for details.

h=interp2(xx,yy,hh,p(:,1),p(:,2),'*linear');

repo: ifem


text:
function [C] = slerp(A,B,t)
  % SLERP Spherically interpolate between to unit vectors. If not unit then
  % this will normalize, slerp and then lerp magnitudes.
  %
  % [C] = slerp(A,B,t)
  %
  % Inputs:
  %   A  #A by 3 list of unit vectors
  %   B  #A by 3 list of unit vectors
  %   t  #t list of values between 0 and 1
  % Outputs:
  %   C  #A by 3 list of interpolated unit vectors
  %

  assert(size(A,1)==size(B,1),'#A != #B');
  if numel(t) == 1
    t = repmat(t,size(A,1),1);
  end
  assert(size(A,1)==numel(t),'#A != #T');

  % Original magnitudes
  Amag = sqrt(sum(A.^2,2));
  Bmag = sqrt(sum(B.^2,2));
  % normalize a and b and keep 0's
  A = bsxfun(@rdivide,A,Amag);
  B = bsxfun(@rdivide,B,Bmag);
  Omega = acos( sum(A.*B,2) );
  Omega(abs(Omega-pi)<1e-8) = pi-10*1e-8;
  sinOmega = sin(Omega);
  C = bsxfun(@times,A,sin((1-t).*Omega)./sinOmega) + ...
      bsxfun(@times,B,sin(t.*Omega)./sinOmega);
  % coincident
  co = abs(Omega)<eps;
  % Converge to typical lerp
  if any(co)
    C(co,:) = A(co,:) + bsxfun(@times,t(co),(B(co,:)-A(co,:)));
  end

  % lerp original magnitudes
  C = bsxfun(@times,C,Amag + t.*(Bmag-Amag));

end

author: alecjacobson

repo: gptoolbox


text:
function new_la = find_good_la(rows,cols,obs,goal,start,old_la,curveset)
%
%find a good look-ahead value to use
%
limits = 0.05;%set threshold
w = extend_world(rows,cols,obs,goal,old_la);
p = calc_pot_value(w);
stpot = p(start.r,start.c,1);
ti = 1;
for i = 1:16
    if curveset.cy(i) <= stpot
        ti = i;
    else
        break;
    end
end
la1 = curveset.cx(ti);
la0 = old_la;
la2 = 500;
for i = ti+1:16
    sig = std(curveset.pty((1:25)+25*(i-1)));
    mu = curveset.cy(i);
    z = (stpot - mu)/sig;
    zVal = normcdf(z);% z table value
    if (zVal <= limits)
        la2 = curveset.cx(i);
    else
        break;
    end
end
a = la2 - la0;
b = la2 - la1;
new_la = la0 + max([a b]);
end


text:
 function axis_equal
%function axis_equal
%|
%| my version of "axis equal"

set(gca, 'DataAspectRatio', [1 1 1])

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
I = imread('tire.tif');
figure(1); imshow(I); figure(2); imhist(I); figure(3); imshow(255-I); figure(4); imhist(255-I);
figure(5); imshow(I); figure(6); imhist(I); figure(7); histeq(I); figure(8); imhist(histeq(I));

author: UtkarshPathrabe


text:
function pvand_test ( )

%*****************************************************************************80
%
%% PVAND_TEST tests PVAND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 November 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PVAND_TEST:\n' );
  fprintf ( 1, '  Solve a Vandermonde linear system A*x=b\n' );

  n = 5;

  for test = 1 : 2

    if ( test == 1 )
      alpha = [ 0.0; 1.0; 2.0; 3.0; 4.0 ];
    elseif ( test == 2 )
      alpha = rand ( n, 1 );
    end

    r8vec_print ( n, alpha, '  Vandermonde vector ALPHA:' );

    a = vand1 ( n, alpha );

    x = [ 5.0; 3.0; 4.0; 1.0; 2.0 ];
    b = a * x;
    r8vec_print ( n, b, '  Right hand side B:' );

    x = pvand ( n, alpha, b );
    r8vec_print ( n, x, '  Solution X:' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function gscatter3(x,y,z,group,clr,sym,siz,doleg,xnam,ynam,znam)
%GSCATTER3  3D Scatter plot with grouping variable
%   gscatter3(x,y,z,group,clr,sym,siz,doleg,xnam,ynam,znam)
%   Designed to work in the exactly same fashion as statistics toolbox's gscatter
%   This function does not require the statistics toolbox installed.
%
%
%   See also GSCATTER, GSCATTER3B
%
%   Copyright 2017 Gustavo Ferraz Trinade.

% Set number of groups
cgroups = unique(group);

cmap = lines(size(cgroups,1));

% Input variables
if (nargin < 5),  clr = lines(max(size(cgroups))); end
if (nargin < 6) || isempty(sym), sym = 'odphs><v^odphs><v^odphs><v^odphs><v^'; end
if (nargin < 7),  siz = 100;           end
if (nargin < 8),  doleg = 'on';        end
if (nargin < 9),  xnam = inputname(1); end
if (nargin < 10), ynam = inputname(2); end
if (nargin < 11), znam = inputname(3); end

% Get current axes
a = gca;
hold(a,'on')

% Plot
for i=1:max(size(cgroups))

    if iscell(cgroups) || ischar(cgroups)
        gi = find(strcmp(group,cgroups(i)));
    else
        gi = find(group == cgroups(i));
    end

    scatter3(a,x(gi),y(gi),z(gi),siz,clr(i,:),'filled',sym(i));
end

% Axes labels and legend (this bit slows down the function)
xlabel(a,xnam);
ylabel(a,ynam);
zlabel(a,znam);

if strcmp(doleg,'on')
    if iscell(cgroups) || ischar(cgroups)
        legend(cgroups')
    else
        legend(num2str(cgroups'))
    end
end

end

author: opencobra

repo: cobratoolbox


text:
function C = mrdivide(A,B)
%MRDIVIDE     Long division A / B
%

% written  12/30/98     S.M. Rump
% modified 10/22/99     S.M. Rump  improvement of error if 1/B exact in few digits
% modfied  02/09/00     S.M. Rump  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
% modified 01/29/10     S.M. Rump  midpoint corrected, thanks to Nozomu Matsuda and Nobito Yamamoto
% modified 06/09/10     S.M. Rump  division w/o error term, thanks to Nozomu Matsuda
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  INTLAB_LONG_BETA = getappdata(0,'INTLAB_LONG_BETA');
  INTLAB_LONG_LOGBETA = getappdata(0,'INTLAB_LONG_LOGBETA');
  INTLAB_LONG_ERROR = getappdata(0,'INTLAB_LONG_ERROR');
  INTLAB_LONG_PRECISION = getappdata(0,'INTLAB_LONG_PRECISION');

  A = long(A);

  if isa(B,'double') & all( B==round(B) ) & all( abs(B)<INTLAB_LONG_BETA )

    % check zero denominator
    if any( B==0 )
      error('long division by zero')
    end

    % Division by integer
    C = A;
    C.sign = C.sign .* sign(B);
    B = abs(B);

    % Compute mantissa
    precA = size(A.mantissa,2);
    C.mantissa(:,1) = floor( A.mantissa(:,1)./B );
    for i=2:precA
      A.mantissa(:,i) = A.mantissa(:,i) + ...
         ( A.mantissa(:,i-1) - B.*C.mantissa(:,i-1) )*INTLAB_LONG_BETA ;
      C.mantissa(:,i) = floor( A.mantissa(:,i)./B );
    end

    % If remainder nonzero, continue division
    rem = A.mantissa(:,precA) - B.*C.mantissa(:,precA);
    if any( rem )
      C.mantissa = ...
         [ C.mantissa zeros(size(C.mantissa,1),INTLAB_LONG_PRECISION-precA+1) ];
      i = precA;
      while any( rem ) & ( i<INTLAB_LONG_PRECISION+1 )
        i = i+1;
        num = rem*INTLAB_LONG_BETA;
        C.mantissa(:,i) = floor( num./B );
        rem = num - B.*C.mantissa(:,i);
      end
    end

    % Omit trailing zeros
    if i<INTLAB_LONG_PRECISION
      C.mantissa = C.mantissa(:,1:i);
    end

    % Compute error
    if INTLAB_LONG_ERROR
      precC = size(C.mantissa,2);
      C.error = errorupdate( -B , C.error , 0 , 1 , rem~=0 , C.exponent-precC );
    end

    % Omit leading zeros
    C = omitleadingzeros(C);

    % Normalize result
    C = normalize(C);
    if INTLAB_LONG_ERROR
      C.error = errornormalize(C.error);
    end

  else   % long division

    % denominator with one extra digit precision
    INTLAB_LONG_PRECISION = INTLAB_LONG_PRECISION + 1;
    B = long(B);
    if INTLAB_LONG_ERROR
      Bmid = mid(B);
      % denominator offset rad(B)^2/mid(B) in double
      N = long2intval(rad(B)).^2./long2intval(Bmid);
      C = addlongerror(Bmid - mid(N),rad(N));     % to be inverted
    else
      Bmid = B;
      C = Bmid;
    end
    precC = size(C.mantissa,2);

    % check zero denominators
    if any( C.mantissa(:,1)==0 )
      error('long division by zero')
    end

    % approximation of inverse for Newton iteration, approx. 52 bits accuracy
    Cinv = C;
    Cinv.exponent = 0;
    Cinv = long(1./long2dble(Cinv));    % inverse of mid(C)
    Cinv.exponent = Cinv.exponent - C.exponent;

    % approximate inverse without error term
    err = longinit('errorterm',0);
    longinit('WithoutErrorTerm',0);

    while 1
      Res = 2 - Cinv*C;
      % check Res==1
      if all( Res.mantissa(:,1)==1 ) & all(all( Res.mantissa(:,2:end)==0 ))
        break
      end
      Cinv = Cinv * Res;
    end

    % reset precision and error term, Cinv~1/C
    longinit(err,0);
    INTLAB_LONG_PRECISION = INTLAB_LONG_PRECISION - 1;

    if INTLAB_LONG_ERROR
      % calculate error
      % error by inversion:  Res_/C  with
      %   Res_ < beta^(1-precRes),   1/C >= 1/C_1*beta^(1-C.exponent)
      % even Res < beta^(-INTLAB_LONG_PRECISION) because Newton iteration
      %   stopped with Res==1 in precision INTLAB_LONG_PRECISION+1; this
      %   is important if 1/C is exactly representable in fewer than
      %   INTLAB_LONG_PRECISION+1 beta-digits (thanks to Kurt Zehetleitner)
      precCinv = size(Cinv.mantissa,2);
      Cmant = C.mantissa(:,1);
      if precC>1
        Cmant = Cmant + C.mantissa(:,2)/INTLAB_LONG_BETA;
      end
      Cinv.error = ...
           errorupdate( -Cmant , 1 , 1-C.exponent-INTLAB_LONG_PRECISION );

      % error introduced by B:  B.error/(sqr(B)-sqr(B.error))
      Bmant = B.mantissa(:,1);
      if size(B.mantissa,2)>1
        Bmant = Bmant + B.mantissa(:,2)/INTLAB_LONG_BETA;
      end
      setround(1)
      N = ( B.error.mant .* INTLAB_LONG_BETA.^(B.error.exp-B.exponent) ) .^ 2;
      setround(-1)
      N = ( Bmant/INTLAB_LONG_BETA ) .^ 2 - N;
      % true denominator  N*beta^(2*B.exponent)

      % check zero denominator (including error)
      if any( N<=0 )
        error('long division by zero')
      end

      % first part:   B.error/(sqr(B)-sqr(B.error))
      Cinv.error = errorupdate( 1 , Cinv.error , 0 , ...
                               -N , B.error , -2*B.exponent );
      Cinv.error = errornormalize(Cinv.error);
    end

    C = A*Cinv;

  end

  setround(rndold)

repo: OpenMAS


text:
function this = orderfields(this, vals);

    for index = 1:length(this)
        this(index).EEG = orderfields(this(index).EEG, vals);
    end;

author: PatternRecognition

repo: OpenBMI


text:
function y = vl_gaussian(x)
% VL_GAUSSIAN  Standard Gaussian density function
%   Y=VL_GAUSSIAN(X) computes the standard (zero mean, unit variance)
%   Gaussian density.
%
%   To obtain the Gaussian density of standard deviation S do
%
%     Y = 1/S * VL_GAUSSIAN(X/S).
%
%   See also: VL_DGAUSSIAN(), VL_DDGAUSSIAN(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

y = 1/sqrt(2*pi)*exp(-0.5*x.^2) ;

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function DoPlotVObj(handles)

%--------Display Parameters
fieldname=fieldnames(handles.Attrh2);
for i=1:length(fieldname)/2
    try
        eval(['SD.' fieldname{i*2} '=' get(handles.Attrh2.(fieldname{i*2}),'String') ';']);
    catch me
        TAttr=get(handles.Attrh2.(fieldname{i*2}),'String');
        eval(['SD.' fieldname{i*2} '=''' TAttr{get(handles.Attrh2.(fieldname{i*2}),'Value')}  ''';']);
    end
end
handles.SD=SD;
%--------End

DoWriteXML2m(DoParseXML(handles.VObjXMLFile),[handles.VObjXMLFile(1:end-3) 'm']);
clear functions;  % remove the M-functions from the memory
[pathstr,name,ext]=fileparts(handles.VObjXMLFile);
eval(['[Obj, VObj]=' name '(0);']);

axes(handles.VObj_axes);
[az,el] = view(gca);
cla(handles.VObj_axes);
for i=1:length(Obj(1,:))
    p(i)=patch(Obj{1,i});
    if isempty(str2num(Obj{2,i}))
        set(p(i),'FaceColor',Obj{2,i},'FaceAlpha',Obj{3,i});
    else
        set(p(i),'FaceColor',str2num(Obj{2,i}),'FaceAlpha',Obj{3,i});
    end
end
view([az,el]);
box(SD.Box);
eval(['grid ' SD.Grid]);
cameratoolbar(SD.CameraTool);
axis image;
set(gca,'XLim',[(-(VObj.XDim-1)/2)*VObj.XDimRes ((VObj.XDim-1)/2)*VObj.XDimRes],...
        'YLim',[(-(VObj.YDim-1)/2)*VObj.YDimRes ((VObj.YDim-1)/2)*VObj.YDimRes],...
        'ZLim',[(-(VObj.ZDim-1)/2)*VObj.ZDimRes ((VObj.ZDim-1)/2)*VObj.ZDimRes]);
set(gca,'YDir','reverse','ZDir','reverse');
xlabel('X');
ylabel('Y');
zlabel('Z');
lighting phong;

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function P = vecRt_r()
% vec(R') = P * r
P = sparse([1,2,3,4,5,6,7,8,9],...
            [1,4,7,2,5,8,3,6,9],...
            [1,1,1,1,1,1,1,1,1],...
            9,9);
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function M_hat = run_mc(params)
  M = params.M;
  A = inexact_alm_mc(M,1e-4,100);
  M_hat = A.U*A.V';
end

author: andrewssobral

repo: lrslibrary


text:
% WAVELET_LAYER_1D Compute the one-dimensional wavelet transform from
% the modulus wavelet coefficients of the previous layer.
%
% Usages
%    [U_phi , U_psi] = wavelet_layer_1d(U, filters)
%
%    [U_phi , U_psi] = wavelet_layer_1d(U, filters, scat_opt)
%
%    [U_phi , U_psi] = wavelet_layer_1d(U, filters, scat_opt, wavelet)
%
% Input
%    U (struct): The input layer to be transformed.
%    filters (cell): The filters of the wavelet transform.
%    scat_opt (struct): The options that are transferred to the wavelet
%    function handle.
%    wavelet (function handle): the wavelet transform function (default is
%    @wavelet_1d).
%
% Output
%    U_phi The coefficients of in, lowpass-filtered (scattering
%       coefficients).
%    U_psi: The wavelet transform coefficients.
%
% Description
%    This function has a pivotal role between WAVELET_1D (which computes a
%    single wavelet transform), and WAVELET_FACTORY_1D (which creates the
%    whole cascade). Given inputs modulus wavelet coefficients
%    corresponding to a layer, WAVELET_LAYER_1D computes the wavelet
%    transform coefficients of the next layer using WAVELET_1D.
% See also
%   WAVELET_1D, WAVELET_FACTORY_1D, WAVELET_LAYER_2D

function [U_phi, U_psi] = wavelet_layer_1dave(U, filters, scat_opt, wavelet)
	if nargin < 3
		scat_opt = struct();
	end

	if nargin < 4
		wavelet = @wavelet_1d;
	end

	calc_U = (nargout>=2);

	[psi_xi,psi_bw,phi_bw] = filter_freq(filters.meta);

	if ~isfield(U.meta, 'bandwidth'), U.meta.bandwidth = 2*pi; end
	if ~isfield(U.meta, 'resolution'), U.meta.resolution = 0; end

	U_phi.signal = {};
	U_phi.meta.bandwidth = [];
	U_phi.meta.resolution = [];
	U_phi.meta.j = zeros(size(U.meta.j,1),0);

	U_psi.signal = {};
	U_psi.meta.bandwidth = [];
	U_psi.meta.resolution = [];
	U_psi.meta.j = zeros(size(U.meta.j,1)+1,0);

	r = 1;
	for p1 = 1:length(U.signal)
		psi_mask = calc_U&(U.meta.bandwidth(p1)>psi_xi);

		scat_opt.x_resolution = U.meta.resolution(p1);
		scat_opt.psi_mask = psi_mask;
		[x_phi, x_psi, meta_phi, meta_psi] = ...
			wavelet(U.signal{p1}, filters, scat_opt);

		U_phi.signal{1,p1} = mean(U.signal{p1}(:));
		U_phi.meta = map_meta(U.meta,p1,U_phi.meta,p1);
		U_phi.meta.bandwidth(1,p1) = meta_phi.bandwidth;
		U_phi.meta.resolution(1,p1) = meta_phi.resolution;

		ind = r:r+sum(psi_mask)-1;
		U_psi.signal(1,ind) = x_psi(1,psi_mask);
		U_psi.meta = map_meta(U.meta,p1,U_psi.meta,ind,{'j'});
		U_psi.meta.bandwidth(1,ind) = meta_psi.bandwidth(1,psi_mask);
		U_psi.meta.resolution(1,ind) = meta_psi.resolution(1,psi_mask);
		U_psi.meta.j(:,ind) = [U.meta.j(:,p1)*ones(1,length(ind)); ...
			meta_psi.j(1,psi_mask)];

		r = r+length(ind);
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [vw, ok] = loadParameterMap(vw, mapPath)
% function [vw, ok] = loadParameterMap(vw, [mapPath])
%
% load and set parameter map in a user specified file
%
% mapPath: optional path to parameter map file. Default: prompt user
%
% rmk, 1/12/99
% djh, 2/21/2001, modified to use dataDir and getPathStrDialog
% ras, 02/04, creates an estimated co, amp, and ph fields if they don't
% exist, b/c mrLoadRet is snotty about using parameter maps alone.
%
% ras, 03/11/04: In working with contrast maps, I've found it useful to have
%  independent statistical threshold and effect size maps -- to restrict
% the map overlay according to the threshold, but to color-code according
% to the effect size. This is analgous to thresholding by coherence, but
% color-coding by amplitude. In this way, for instance, one could leave the
% mapWin at full range, slide the cothresh sliand see both highly correlated
% and anticorrelated regions for a given contrast.
% So, I've made some param map files that have a 'co' field saved with the
% map and mapName vars. This field, when loaded, would replace the co field
% for the corAnal.
% ARW: It didn't seem right that we estimated co (or had it saved in the
% param file). Instead, we fill co with ones so that it behaves correctly
% ras: Yeah, it's a hack, b/c the variable doesn't generally represent coherence,
% but for event-related analyses, it's often pretty useful to have the
% two different ways of thresholding. While a long-term solution is to have
% a more generalized display method (which is in the nascent mrVista2
% viewer), it's still nice to have it for now.
%
% ras, 11/15/05: added ability to load map parameters, such as the color
% map and clip modes, which may be saved along with the map.
%
% rkl, 04/17/17: Fixed a bug: the UI display mode did not get set to 'map'
% when loading a nifti as a parameter map (though everything works smoothly
% when loading a .mat file as parameter map). This resulted in the slider not
% displaying the full range of values.
if ~exist('mapPath','var')
    mapPath = getPathStrDialog(dataDir(vw),'Choose parameter map file name',{'*.mat; *.nii; *.nii.gz'});
    if isempty(mapPath), ok = false; return; end
end

% if the map path doesn't exist, but a map file exists in the view's
% data dir with that name, modify the path accordingly:
if ~check4File(mapPath)
    if check4File( fullfile(dataDir(vw), mapPath) )
        mapPath = fullfile(dataDir(vw), mapPath);
    end
end

[~,name,extension] = fileparts(mapPath);

if strcmp(extension, '.mat')
    % check4File only works with .mat files:
    if ~check4File(mapPath)
        warning(['No ',mapPath,' file']) %#ok<WNTAG>
        ok = 0;
        return
    end

else
    % We still need to check whether the file exists:
    if ~exist(mapPath, 'file')
        warning(['No ',mapPath,' file']) %#ok<WNTAG>
        ok = 0;
        return
    end
end

verbose = prefsVerboseCheck;
if verbose
    fprintf('[%s]: Loading Parameter Map: %s \n', mfilename, mapPath);
end

% If it's a nifti file:
if strcmp(extension, '.nii') || strcmp(extension,'.gz')

    % Read in the nifti
    nii = niftiRead(mapPath);

    % Apply the canonical xform for nifti data
    switch lower(viewGet(vw, 'viewType'))
        case 'inplane'
            nii = niftiApplyAndCreateXform(nii,'Inplane');

            % Set the data into the view structure, converting a 4D array
            % into a cell array of 3D arrays
            niData = niftiGet(nii,'Data');
            n = size(niData,4);
            mapData = cell(1,n);
            for ii = 1:n
                mapData{ii} = niData(:,:,:,ii);
            end
            vw = viewSet(vw, 'map', mapData);

            % To deal with with the double-extension in .nii.gz files
            dot_idx = find(name=='.', 1);
            if ~isempty(dot_idx), name = name(1:find(name=='.')-1); end

            % It will be named according to the name of the nifti file:
            vw = viewSet(vw, 'map Name', name);

            % Save it in the right place for the DATATYPE:
            saveParameterMap(vw)

            mapPath = fullfile(dataDir(vw),[vw.mapName,'.mat']);

            load(mapPath);

            vw = setParameterMap(vw,map,mapName);

        case {'gray' 'volume'}
            vw = nifti2functionals(vw, mapPath);

            % refresh (if this isn't a hidden view -- no ui struct)
            if isfield(vw, 'ui')
                vw = setDisplayMode(vw,'map');
            end

            return;
        otherwise

    end

end

% Load it
load(mapPath);

% Set the parameters
vw = setParameterMap(vw,map,mapName);

nScans = viewGet(vw, 'numScans');
if exist('co','var')
    for scan = 1:nScans
        if length(co) >= scan && ~isempty(co{scan}) %#ok<USENS>
            vw.co{scan} = co{scan};
        end
    end
end

% are the units for the map saved?
% (due to inconsistent coding, there are 2 possible variable names.
%  'units' or 'mapUnits'. 'mapUnits' is preferred, but for
%  back-compatibility, I leave in a check for 'units' as well.)
if exist('mapUnits', 'var')
    vw = viewSet(vw, 'mapUnits', mapUnits);
elseif exist('units', 'var')
    vw = viewSet(vw, 'mapUnits', units);
else
        vw = viewSet(vw, 'mapUnits', ''); % reset to empty
end

% if parameters about the view mode are saved, load these in as well
if exist('cmap','var'),     vw.ui.mapMode.cmap = cmap; end
if exist('clipMode','var'), vw.ui.mapMode.clipMode = clipMode; end
if exist('numColors','var'),vw.ui.mapMode.numColors = numColors; end
if exist('numGrays','var'), vw.ui.mapMode.numGrays = numGrays; end

ok = 1;

% refresh (if this isn't a hidden view -- no ui struct)
if isfield(vw, 'ui')
    vw = setDisplayMode(vw,'map');
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%--- help for tabular/summary ---
%
% SUMMARY Print summary of a table or a timetable.
%    SUMMARY(T) prints a summary of T and the variables that it contains. If
%    T is a table, then SUMMARY displays the description from
%    T.Properties.Description followed by a summary of the table variables.
%    If T is a timetable, then SUMMARY additionally displays a summary of
%    the row times.
%
%    S = SUMMARY(T) returns a summary of T as a structure. Each field of S
%    contains a summary for the corresponding variable in T. If T is a
%    timetable, S contains an additional field for a summary of the row
%    times.
%
%    Other functions named summary
%
%       categorical/summary      dataset/summary        mcmc/summary
%       codistributed/summary    distributed/summary    tall/summary
%

author: jmaih


text:
function y = ainv_x ( x, R, k )

%*****************************************************************************80
%
%% AINV_X computes inverse(A)*x.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(*), the vector to be multiplied by inverse(A).
%
%    Input, character R, the first argument to the NUMGRID command.
%
%    Input, integer K, the order of the grid.
%
%    Output, real Y, the value of inverse(A)*X.
%
  y = ( delsq ( numgrid ( R, k ) ) ) \ x;

  return
end


text:
function [ a, b ] = p02_lim ( dim_num )

%*****************************************************************************80
%
%% P02_LIM returns the integration limits for problem 02.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Output, real A(DIM_NUM), B(DIM_NUM), the lower and upper
%    limits of integration.
%
  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  return
end


text:
function rule_num = triangle_ncc_rule_num ( )

%*****************************************************************************80
%
%% TRIANGLE_NCC_RULE_NUM returns the number of NCC rules available.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Peter Silvester,
%    Symmetric Quadrature Formulae for Simplexes,
%    Mathematics of Computation,
%    Volume 24, Number 109, January 1970, pages 95-100.
%
%  Parameters:
%
%    Output, integer RULE_NUM, the number of rules available.
%
  rule_num = 9;

  return
end


text:
function rtk=udclk_ppp(rtk)

global glc
opt=rtk.opt; VAR_CLK=60^2; CLIGHT=glc.CLIGHT; tt=abs(rtk.tt);
navsys=opt.navsys; mask=rtk.mask; prn=0;

% for single-system except GPS
if strcmp(navsys,'R')
    dtr=rtk.sol.dtr(2);
    if abs(dtr)<=1e-16,dtr=1e-16;end
    rtk=initx(rtk,CLIGHT*dtr,VAR_CLK,rtk.ic+2);

    % for GLONASS icb
    if opt.gloicb==glc.GLOICB_LNF
        if rtk.x(rtk.iicb+1)==0
            rtk=initx(rtk,0.1,VAR_CLK,rtk.iicb+1);
        end
    elseif opt.gloicb==glc.GLOICB_QUAD
        if rtk.x(rtk.iicb+1)==0
            rtk=initx(rtk,0.1,VAR_CLK,rtk.iicb+1);
        end
        if rtk.x(rtk.iicb+2)==0
            rtk=initx(rtk,0.1,VAR_CLK,rtk.iicb+2);
        end
    end
    return
elseif strcmp(navsys,'E')
    dtr=rtk.sol.dtr(3);
    if abs(dtr)<=1e-16,dtr=1e-16;end
    rtk=initx(rtk,CLIGHT*dtr,VAR_CLK,rtk.ic+3);
    return
elseif strcmp(navsys,'C')
    dtr=rtk.sol.dtr(4);
    if abs(dtr)<=1e-16,dtr=1e-16;end
    rtk=initx(rtk,CLIGHT*dtr,VAR_CLK,rtk.ic+4);
    return
elseif strcmp(navsys,'J')
    dtr=rtk.sol.dtr(5);
    if abs(dtr)<=1e-16,dtr=1e-16;end
    rtk=initx(rtk,CLIGHT*dtr,VAR_CLK,rtk.ic+5);
    return
end

% for GPS or multi-system
dtr=rtk.sol.dtr(1);
if abs(dtr)<=1e-16,dtr=1e-16;end
rtk=initx(rtk,CLIGHT*dtr,VAR_CLK,rtk.ic+1);
for i=2:glc.NSYS
    if mask(i)==0,continue;end
    if i==glc.SYS_GLO
        dtr=rtk.sol.dtr(2);
        if rtk.x(rtk.ic+2)==0
            if abs(dtr)<=1e-16,dtr=1e-16;end
            rtk=initx(rtk,CLIGHT*dtr,VAR_CLK,rtk.ic+2);
        else
            rtk.P(rtk.ic+i,rtk.ic+i)=rtk.P(rtk.ic+i,rtk.ic+i)+prn^2*tt;
        end
        % for GLONASS icb
        if opt.gloicb==glc.GLOICB_LNF
            if rtk.x(rtk.iicb+1)==0
                rtk=initx(rtk,0.1,VAR_CLK,rtk.iicb+1);
            end
        elseif opt.gloicb==glc.GLOICB_QUAD
            if rtk.x(rtk.iicb+1)==0
                rtk=initx(rtk,0.1,VAR_CLK,rtk.iicb+1);
            end
            if rtk.x(rtk.iicb+2)==0
                rtk=initx(rtk,0.1,VAR_CLK,rtk.iicb+2);
            end
        end
    else
        dtr=rtk.sol.dtr(i);
        if rtk.x(rtk.ic+i)==0
            if abs(dtr)<=1e-16,dtr=1e-16;end
            rtk=initx(rtk,CLIGHT*dtr,VAR_CLK,rtk.ic+i);
        else
            rtk.P(rtk.ic+i,rtk.ic+i)=rtk.P(rtk.ic+i,rtk.ic+i)+prn^2*tt;
        end
    end
end

return

repo: GINav


text:
classdef GlobalBaseInertialFrame < AbstractReferenceFrame
    %GlobalBaseInertialFrame Summary of this class goes here
    %   Detailed explanation goes here

    properties
        celBodyData
    end

    properties(Constant)
        typeEnum = ReferenceFrameEnum.BodyCenteredInertial
    end

    methods
        function obj = BodyCenteredInertialFrame(celBodyData)
            obj.celBodyData = celBodyData;
        end

        function [posOffsetOrigin, velOffsetOrigin, angVelWrtOrigin, rotMatToInertial] = getOffsetsWrtInertialOrigin(obj, time, ~, ~)
            posOffsetOrigin = repmat([0;0;0], [1, length(time)]);
            velOffsetOrigin = repmat([0;0;0], [1, length(time)]);
            [angVelWrtOrigin, rotMatToInertial] = obj.getAngVelWrtOriginAndRotMatToInertial(time, [], []);
        end

        function rotMatToInertial = getRotMatToInertialAtTime(~, time, ~, ~)
            rotMatToInertial = repmat(eye(3), [1, 1, length(time)]);
        end

        function [angVelWrtOrigin, rotMatToInertial] = getAngVelWrtOriginAndRotMatToInertial(obj, time, vehElemSet, bodyInfoInertialOrigin)
            angVelWrtOrigin = repmat([0;0;0], [1, length(time)]);
            rotMatToInertial = obj.getRotMatToInertialAtTime(time, vehElemSet, bodyInfoInertialOrigin);
        end

        function tf = frameOriginIsACelBody(obj)
            tf = false;
        end

        function bodyInfo = getOriginBody(~)
            error('Not supported.');
        end

        function setOriginBody(~, ~)
            %nothing to do
        end

        function nameStr = getNameStr(~)
            nameStr = sprintf('Global Base Inertial Frame');
        end

        function newFrame = editFrameDialogUI(obj, ~)
            newFrame = obj;
        end
    end

    methods
        function bool = eq(A,B)
            bool = eq@handle(A,B);
        end

        function bool = ne(A,B)
            bool = ne@handle(A,B);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function a = circulant ( m, n, x )

%*****************************************************************************80
%
%% CIRCULANT returns the CIRCULANT matrix.
%
%  Formula:
%
%    K = 1 + mod ( J-I, N )
%    A(I,J) = X(K)
%
%  Example:
%
%    M = 4, N = 4, X = ( 1, 2, 3, 4 )
%
%    1  2  3  4
%    4  1  2  3
%    3  4  1  2
%    2  3  4  1
%
%    M = 4, N = 5, X = ( 1, 2, 3, 4, 5 )
%
%    1  2  3  4  5
%    5  1  2  3  4
%    4  5  1  2  3
%    3  4  5  1  2
%
%    M = 5, N = 4, X = ( 1, 2, 3, 4 )
%
%    1  2  3  4
%    5  1  2  3
%    4  5  1  2
%    3  4  5  1
%    1  2  3  4
%
%  Discussion:
%
%    Westlake lists the following "special" circulants:
%
%      B2, X = ( T^2, 1, 2, ..., T, T+1, T, T-1, ..., 1 ),
%      with T = ( N - 2 ) / 2;
%
%      B3, X = ( N+1, 1, 1, ..., 1 );
%
%      B5, X = ( 1, 2, 3, ..., N ).
%
%  Rectangular Properties:
%
%    The product of two circulant matrices is a circulant matrix.
%
%    The transpose of a circulant matrix is a circulant matrix.
%
%    A circulant matrix C, whose first row is (c1, c2, ..., cn), can be
%    written as a polynomial in the upshift matrix U:
%
%      C = c1 * I + c2 * U + c3 * U**2 + ... + cn * U**n-1.
%
%    A is a circulant: each row is shifted once to get the next row.
%
%  Square Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is Toeplitz: constant along diagonals.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    A commutes with any other circulant matrix.
%
%    A is normal.
%
%    The transpose of A is also a circulant matrix.
%
%    The inverse of A is also a circulant matrix.
%
%    The Fourier matrix is the eigenvector matrix for every circulant matrix.
%
%    Because the Fourier matrix F diagonalizes A, the inverse (or
%    pseudoinverse, if any LAMBDA is zero) can be written
%
%      inverse ( A ) = (F*) * 1/LAMBDA * F
%
%    A is symmetric if, for all I, X(I+1) = X(N-I+1).
%
%    If R is an N-th root of unity, that is, R is a complex number such
%    that R**N = 1, then
%
%      Y = X(1) + X(2)*R + X(3)*R**2 + ... + X(N)*R**(N-1)
%
%    is an eigenvalue of A, with eigenvector
%
%     ( 1, R, R**2, ..., R**(N-1) )
%
%    and left eigenvector
%
%      ( R**(N-1), R**(N-2), ..., R**2, R, 1 ).
%
%    Although there are exactly N distinct roots of unity, the circulant
%    may have repeated eigenvalues, because of the behavior of the polynomial.
%    However, the matrix is guaranteed to have N linearly independent
%    eigenvectors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Davis,
%    Circulant Matrices,
%    John Wiley, 1979, QA188.D37.
%
%    Robert Gregory, David Karney,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, New York, 1969, page 22,
%    LC: QA263.G68.
%
%    Joan Westlake,
%    Test Matrix A24,
%    A Handbook of Numerical Matrix Inversion and Solution of Linear Equations,
%    John Wiley, 1968.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of A.
%
%    Input, real X(N), the values in the first row of A.
%
%    Output, real A(M,N), the matrix.
%
  for i = 1 : m
    for j = 1 : n

      k = 1 + i4_modp ( j-i, n );
      a(i,j) = x(k);

    end
  end

  return
end


text:
function pwl_interp_2d_scattered_test ( )

%*****************************************************************************80
%
%% PWL_INTERP_2D_SCATTERED_TEST tests the PWL_INTERP_2D_SCATTERED library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 October 2012
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../pwl_interp_2d_scattered' )
  addpath ( '../test_interp_2d' )
  addpath ( '../r8lib' )

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PWL_INTERP_2D_SCATTERED_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test PWL_INTERP_2D_SCATTERED.\n' );
  fprintf ( 1, '  The R8LIB library is needed.\n' );
  fprintf ( 1, '  This test also needs the TEST_INTERP_2D library.\n' );

  pwl_interp_2d_scattered_test01 ( );
  pwl_interp_2d_scattered_test02 ( );
%
%  Numerical tests.
%
  prob_num = f00_num ( );

  for prob = 1 : prob_num
    pwl_interp_2d_scattered_test03 ( prob );
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PWL_INTERP_2D_SCATTERED_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../pwl_interp_2d_scattered' )
  rmpath ( '../test_interp_2d' )
  rmpath ( '../r8lib' )

  return
end


function plot3Dmodel(MatFileName)
%PLOT3DMODEL(MatFileName)
% Plots the given 3d model, test for animation function
%

load(MatFileName)

%% Check the results
% Get maximum dimension to plot the circles afterwards
AC_DIMENSION = max(max(sqrt(sum(Model3D.Aircraft(1).stl_data.vertices.^2,2))));
for i=1:length(Model3D.Control)
    AC_DIMENSION = max(AC_DIMENSION,max(max(sqrt(sum(Model3D.Control(i).stl_data.vertices.^2,2)))));
end
% Define the figure properties
AX = axes('position',[0.0 0.0 1 1]);
axis off
scrsz = get(0,'ScreenSize');
set(gcf,'Position',[scrsz(3)/40 scrsz(4)/12 scrsz(3)/2*1.0 scrsz(3)/2.2*1.0],'Visible','on');
set(AX,'color','none');
axis('equal')
hold on;
cameratoolbar('Show')
% Circles around the aircraft transformation group handles
euler_hgt(1)  = hgtransform('Parent',           AX, 'tag', 'OriginAxes');
euler_hgt(2)  = hgtransform('Parent', euler_hgt(1), 'tag', 'roll_disc');
euler_hgt(3)  = hgtransform('Parent', euler_hgt(1), 'tag', 'pitch_disc');
euler_hgt(4)  = hgtransform('Parent', euler_hgt(1), 'tag', 'heading_disc');
euler_hgt(5)  = hgtransform('Parent', euler_hgt(2), 'tag', 'roll_line');
euler_hgt(6)  = hgtransform('Parent', euler_hgt(3), 'tag', 'pitch_line');
euler_hgt(7)  = hgtransform('Parent', euler_hgt(4), 'tag', 'heading_line');
% Plot objects
% -------------------------------------------------------------------------
% Plot airframe
for i = 1:length(Model3D.Aircraft)
    AV = patch(Model3D.Aircraft(i).stl_data,  'FaceColor',        Model3D.Aircraft(i).color, ...
        'EdgeColor',        'none',        ...
        'FaceLighting',     'gouraud',     ...
        'AmbientStrength',   0.15);
end
CONT(length(Model3D.Control))=0;
% Plot controls
for i=1:length(Model3D.Control)
    CONT(i) = patch(Model3D.Control(i).stl_data,  'FaceColor',        Model3D.Control(i).color, ...
        'EdgeColor',        'none',        ...
        'FaceLighting',     'gouraud',     ...
        'AmbientStrength',  0.15);
    % Plot the rotation point and the rotation axis of each control
    % (double-check correct implementation and rotation direction of each
    % control surface)
    p = Model3D.Control(i).rot_point;
    vect = Model3D.Control(i).rot_vect;
    plot3(p(1)+[0, AC_DIMENSION*vect(1)/2], p(2)+[0, AC_DIMENSION*vect(2)/2], p(3)+[0, AC_DIMENSION*vect(3)/2], 'b-o', 'MarkerSize', 10, 'LineWidth', 2);
end
% Fixing the axes scaling and setting a nice view angle
axis('equal');
axis([-1 1 -1 1 -1 1] * 2.0 * AC_DIMENSION)
set(gcf,'Color',[1 1 1])
axis off
view([30 10])
zoom(2.0);
% Add a camera light, and tone down the specular highlighting
camlight('left');
material('dull');

% --------------------------------------------------------------------
% Define the radius of the sphere
R = 1.0 * AC_DIMENSION;

% Outer circles
phi = (-pi:pi/36:pi)';
D1 = [sin(phi) cos(phi) zeros(size(phi))];
HP(1) = plot3(R*D1(:,1),R*D1(:,2),+R*D1(:,3),'Color','b','tag','Zplane','Parent',euler_hgt(4));
HP(2) = plot3(R*D1(:,2),R*D1(:,3),+R*D1(:,1),'Color',[0 0.8 0],'tag','Yplane','Parent',euler_hgt(3));
HP(3) = plot3(R*D1(:,3),R*D1(:,1),+R*D1(:,2),'Color','r','tag','Xplane','Parent',euler_hgt(2));

% +0,+90,+180,+270 Marks
S = 0.95;
phi = -pi+pi/2:pi/2:pi;
D1 = [sin(phi); cos(phi); zeros(size(phi))];
plot3([S*R*D1(1,:); R*D1(1,:)],[S*R*D1(2,:); R*D1(2,:)],[S*R*D1(3,:); R*D1(3,:)],'Color','b','tag','Zplane','Parent',euler_hgt(4));
plot3([S*R*D1(2,:); R*D1(2,:)],[S*R*D1(3,:); R*D1(3,:)],[S*R*D1(1,:); R*D1(1,:)],'Color',[0 0.8 0],'tag','Yplane','Parent',euler_hgt(3));
plot3([S*R*D1(3,:); R*D1(3,:)],[S*R*D1(1,:); R*D1(1,:)],[S*R*D1(2,:); R*D1(2,:)],'Color','r','tag','Xplane','Parent',euler_hgt(2));
text(R*1.05*D1(1,:),R*1.05*D1(2,:),R*1.05*D1(3,:),{'N','E','S','W'},'Fontsize',9,'color',[0 0 0],'HorizontalAlign','center','VerticalAlign','middle');

% +45,+135,+180,+225,+315 Marks
S = 0.95;
phi = -pi+pi/4:2*pi/4:pi;
D1 = [sin(phi); cos(phi); zeros(size(phi))];
plot3([S*R*D1(1,:); R*D1(1,:)],[S*R*D1(2,:); R*D1(2,:)],[S*R*D1(3,:); R*D1(3,:)],'Color','b','tag','Zplane','Parent',euler_hgt(4));
HT = text(R*1.05*D1(1,:),R*1.05*D1(2,:),R*1.05*D1(3,:),{'NW','NE','SE','SW'},'Fontsize',8,'color',[0 0 0],'HorizontalAlign','center','VerticalAlign','middle');

% 10 deg sub-division marks
S = 0.98;
phi = -[0:10:90 80:-10:0 -10:-10:-90 -80:10:0];
PHI_TEXT{length(phi)}='';
for i=1:length(phi)
    PHI_TEXT{i} = num2str(phi(i));
end
theta_t = -[0:10:90 80:-10:0 -10:-10:-90 -80:10:0];
THETA_TEXT{length(theta_t)}='';
for i=1:length(theta_t)
    THETA_TEXT{i} = num2str(theta_t(i));
end
phi = -180:10:180;
phi = phi*pi/180;
D1 = [sin(phi); cos(phi); zeros(size(phi))];
plot3([S*R*D1(1,:); R*D1(1,:)],[S*R*D1(2,:); R*D1(2,:)],[S*R*D1(3,:); R*D1(3,:)],'Color','b','tag','Zplane','Parent',euler_hgt(4));
plot3([S*R*D1(2,:); R*D1(2,:)],[S*R*D1(3,:); R*D1(3,:)],[S*R*D1(1,:); R*D1(1,:)],'Color',[0 0.8 0],'tag','Yplane','Parent',euler_hgt(3));
plot3([S*R*D1(3,:); R*D1(3,:)],[S*R*D1(1,:); R*D1(1,:)],[S*R*D1(2,:); R*D1(2,:)],'Color','r','tag','Xplane','Parent',euler_hgt(2));

% Plot guide lines
HL(1) = plot3([-R R],[0 0],[0 0],'b-','tag','heading_line','parent',euler_hgt(7));
HL(2) = plot3([-R R],[0 0],[0 0],'g-','tag','pitch_line','parent',euler_hgt(6),'color',[0 0.8 0]);
HL(3) = plot3([0 0],[-R R],[0 0],'r-','tag','roll_line','parent',euler_hgt(5));
end


text:
function average_velocity(nlp, bounds)
    % constraints for impact velocities

    domain = nlp.Plant;
    x = domain.States.x;

    % average step velocity
    v_lb = bounds.constrBounds.averageVelocity.lb;
    v_ub = bounds.constrBounds.averageVelocity.ub;
    x0 = x;
    xf = SymVariable('xf',size(x));
    T  = SymVariable('t',[2,1]);
    v_avg = [(xf(1)-x0(1))./(T(2)-T(1))
        (xf(2)-x0(2))./(T(2)-T(1))];
    v_avg_fun = SymFunction(['avgStepVelocity_',domain.Name],v_avg,{T, x0, xf});
    x0_var = nlp.OptVarTable.x(1);
    xf_var = nlp.OptVarTable.x(end);
    t_var  = nlp.OptVarTable.T(1);
    v_avg_cstr = NlpFunction('Name',v_avg_fun.Name,...
        'Dimension',2,'lb',v_lb, 'ub',v_ub,'Type','Nonlinear',...
        'SymFun',v_avg_fun,'DepVariables',[t_var;x0_var;xf_var]);
    addConstraint(nlp, 'avgStepVelocity', 'first', v_avg_cstr);

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function X = real(X)
%REAL (overloaded)

X.basis = real(X.basis);
X = clean(X);
if isa(X,'sdpvar')
    X.conicinfo = [0 0];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [ bboxes ] = scale_bboxes( bboxes, scale_factor )
% scale_bboxes: it scales the set bounding boxes bboxes by the scale_factor
% factor.
%
% INPUT:
% 1) bboxes: a N x 4 array with the input bounding box coordinates in the
% form of [x0,y0,x1,y1] (where (x0,y0) is the top-left corner and (x1,y1)
% the bottom left corner)
% 2) scale_factor: a 1 x 1 or 2 x 1 array with the scaling factor of the
% bounding boxes. If scale_factor is a 1 x 1 array then the same scaling
% factor will be applied on both the x and y axis. If scale_factor is a
% 2 x 1 (or 1 x 2) array then the bounding boxes will be scaled across the
% y dimension by scale_factor(1) and across the x dimension by
% scale_factor(2).
%
% OUTPUT:
% 1) bboxes: a N x 4 array with the output bounding box coordinates in the
% form of [x0,y0,x1,y1] (where (x0,y0) is the top-left corner and (x1,y1)
% the bottom left corners.
%
% This file is part of the code that implements the following paper:
% Title      : "LocNet: Improving Localization Accuracy for Object Detection"
% Authors    : Spyros Gidaris, Nikos Komodakis
% Institution: Universite Paris Est, Ecole des Ponts ParisTech
% ArXiv link : http://arxiv.org/abs/1511.07763
% code       : https://github.com/gidariss/LocNet
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2016 Spyros Gidaris
%
% Title     : "LocNet: Improving Localization Accuracy for Object Detection"
% ArXiv link: http://arxiv.org/abs/1511.07763
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

if numel(scale_factor) == 1, scale_factor(2) = scale_factor(1); end
assert(numel(scale_factor) == 2);
scale_factor = single(scale_factor);

bboxes_center      = [(bboxes(:,1)+bboxes(:,3)), (bboxes(:,2)+bboxes(:,4))]/2;
bboxes_width_half  = (bboxes(:,3) - bboxes(:,1))/2;
bboxes_width_half  = bboxes_width_half * scale_factor(2);

bboxes_height_half = (bboxes(:,4) - bboxes(:,2))/2;
bboxes_height_half = bboxes_height_half * scale_factor(1);

bboxes = round([bboxes_center(:,1) - bboxes_width_half, ...
                bboxes_center(:,2) - bboxes_height_half, ...
                bboxes_center(:,1) + bboxes_width_half, ...
                bboxes_center(:,2) + bboxes_height_half]);
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function [h,X,Y] = ellipseupdate(h,a,b,x0,y0,phi,lineStyle)
%ELLIPSEDRAW can draw an arbitrary ellipse with given parameters.
%   The properties of that ellipse plot can be customized
%   by setting the ellipse handle.
%
%       hEllipse = ellipsedraw(a,b,x0,y0,phi,lineStyle)
%
%   Input parameters:
%       a           Value of the major axis
%       b           Value of the minor axis
%       x0          Abscissa of the center point of the ellipse
%       y0          Ordinate of the center point of the ellipse
%       phi         Angle between x-axis and the major axis
%       lineStyle   Definition of the plotted line style
%
%   Output:
%       hEllipse    Handle of the ellipse
%
%   Simple usage:
%       ellipsedraw(5,3);
%       ellipsedraw(5,3,'g--');
%       ellipsedraw(5,3,pi/4);
%
%   Complete usage:
%       h = ellipsedraw(5,3,1,-2,pi/4,'r-.');
%       set(h,'LineWidth',2);

% Designed by: Lei Wang, <WangLeiBox@hotmail.com>, 25-Mar-2003.
% Last Revision: 01-Apr-2003.
% Dept. Mechanical & Aerospace Engineering, NC State University.
% Copyright (c)2003, Lei Wang <WangLeiBox@hotmail.com>
%$Revision: 1.1 $  $ 4/1/2003 5:42:24 PM $

if (nargin < 3)||(nargin > 7),
    error('Please see help for INPUT DATA.');

elseif nargin == 3
    x0 = 0;     y0 = 0;
    phi = 0;    lineStyle = 'b-';

elseif nargin == 4
    if ischar(x0) == 1
        lineStyle = x0;
        x0 = 0; y0 = 0;
        phi = 0;
    else
        phi = x0;
        x0 = 0; y0 = 0;
        lineStyle = 'b-';
    end

elseif nargin == 5
    phi = 0;    lineStyle = 'b-';

elseif nargin == 6
    lineStyle = 'b-';
end

theta = [-0.03:0.01:2*pi];

% Parametric equation of the ellipse
%----------------------------------------
 x = a*cos(theta);
 y = b*sin(theta);

% Coordinate transform
%----------------------------------------
 X = cos(phi)*x - sin(phi)*y;
 Y = sin(phi)*x + cos(phi)*y;
 X = X + x0;
 Y = Y + y0;

% Plot the ellipse
%----------------------------------------
set(h,'xdata',X,'ydata',Y);
%axis equal;

author: kristinbranson

repo: JAABA


text:
% SENSE_recon_demo.m
% - demonstrates use of SENSE_recon.m, which implements the SENSE
% reconstruction
% - figure shows how Tikhonov regularization mitigates the increase in error
% as noise variance increases
% - applies complex AWGN to kspace
% - image must be even both directions
%
% 2012-06-08, Mai Le
% 2013-01-02 Mai Le, tweaks

% comparison = 1 will demonstrate effect of masks and regularization on MSE,
% but demo takes much longer. comparison = 0 runs masked SENSE w/o
% regularization
comparison = 0;

% number of coils
nc = 4;
% factor of undersampling
np = 2;
% dimension to undersample
for reduced_dim = [1 2]
    display(['SENSE reduced in dim ' num2str(reduced_dim)]);

    % set up "true" brain image
    %f.dir = [path_find_dir('mri') '/../data/mri/'];
    %f.xtrue = [f.dir 'brainweb_t1.jpg'];
    %mag_true = double(imread(f.xtrue))'; % true image magnitude
    mag_true = 25*ellipse_im(256);
    % requires image to be multiple of np in dimension of undersampling
    % and even for generating smap
    if (reduced_dim == 1)
        mag_true = mag_true(1:floor(size(mag_true,1)/(np*2))*(np*2),:);
    elseif (reduced_dim == 2)
        mag_true = mag_true(:,1:floor(size(mag_true,2)/(np*2))*(np*2));
    else
        fail 'bug'
    end

    dims = size(mag_true);
    % introduce planar phase
    [xx,yy] = ndgrid(-dims(1)/2:dims(1)/2-1,-dims(2)/2:dims(2)/2-1);
    ph_max = pi/5;
    ph = ph_max*xx/(dims(1)/2)+ph_max*yy/(dims(2)/2);
    image = mag_true.*exp(1i*ph);

    % generate sensitivity maps
    smap = mri_sensemap_sim('nx',dims(1),'ny',dims(2),'ncoil',nc,'rcoil',400);
    close;

    im_rep = repmat(image,[1 1 nc]);

    % apply sensitivity maps
    mapped_im = im_rep.*smap;

    % throw away k-space data
    im_fft = zeros(size(mapped_im));
    if reduced_dim == 1
        reduced_fft = zeros(dims(1)/np,dims(2),nc);
    elseif reduced_dim == 2
        reduced_fft = zeros(dims(1),dims(2)/np,nc);
    else
        fail 'bug'
    end
    for ii = 1:nc
        im_fft(:,:,ii) = fft2(mapped_im(:,:,ii));
        if reduced_dim == 1
            reduced_fft(:,:,ii) = im_fft(1:np:end,:,ii);
        elseif reduced_dim == 2
            reduced_fft(:,:,ii) = im_fft(:,1:np:end,ii);
        else
            fail 'bug'
        end
    end

    %% construct mask
    thresh = 0.005;
    thresh_mask = any(abs(mapped_im)>thresh,3);

    %% perform SENSE reconstruction for varying levels of noise
    f.sigmas = 0:50:550; % complex AWGN sdtd dev
    clear recon
    recon(1).im = zeros(dims(1),dims(2),length(f.sigmas));
    recon(1).err = zeros(length(f.sigmas),1);
    if comparison
        recon = repmat(recon,4,1);
    end
    nrms_fun = @(x,y) 100 * nrms(x(:), y(:));
    for ii = 1:length(f.sigmas)
        % introduce complex AWGN
        sigma = f.sigmas(ii);
        noisy_reduced_fft = reduced_fft + sigma*randn(size(reduced_fft)) + ...
            1i*sigma*randn(size(reduced_fft));

        % SENSE reconstruction
        recon(1).im(:,:,ii) = SENSE_recon(smap, noisy_reduced_fft, ...
            np, 'mask', thresh_mask, 'direction', reduced_dim);
        if comparison
            recon(2).im(:,:,ii) = SENSE_recon(smap, noisy_reduced_fft, ...
                np, 'direction', reduced_dim);
            recon(3).im(:,:,ii) = SENSE_recon(smap, noisy_reduced_fft, ...
                np, 'reg', 'Tikhonov', 'direction', reduced_dim);
            recon(4).im(:,:,ii) = SENSE_recon(smap, noisy_reduced_fft, ...
                np, 'mask', thresh_mask, 'reg', 'Tikhonov', 'direction', reduced_dim);
        end
        for num = 1:(1+3*comparison) % do all 4 if comparison
            recon(num).err(ii) = nrms_fun(image, recon(num).im(:,:,ii));
        end
    end
    %% plot error as a function of noise variance
    figure; plot(f.sigmas.^2, recon(1).err, 'bo');
    if (comparison)
        hold on; plot(f.sigmas.^2, recon(2).err,'ko');
        plot(f.sigmas.^2, recon(3).err, 'rx');
        plot(f.sigmas.^2, recon(4).err, 'gx');
        legend('no reg, mask', 'no reg, no mask', 'tikhonov, no mask', 'tikhonov, mask');
    end
    xlabel('variance of complex AWGN in kspace');
    ylabel('RMS error in reconstructed images');
    if (comparison)
        title('effect of Tikhonov regularization on noisy SENSE reconstruction')
    else
        title('RMSE as function of noise variance');
    end
    %% plot reconstructed image at sigma(ii)
    ii = 4; % change as desired
    clim = minmax(mag_true)';
    if (comparison)
        im(1, real(image));
        im(2, imag(image));
        titlef('imag of original image, noise variance: %d', f.sigmas(ii).^2);
        im(5, real(recon(1).im(:,:,ii)), clim);
        title('real of recon image, no reg');
        im(6, imag(recon(1).im(:,:,ii)), clim);
        title('imag of recon image, no reg');
        im(7, abs(real(image)-real(recon(1).im(:,:,ii)))+ ...
            abs(imag(image)-imag(recon(1).im(:,:,ii))), [-9 9]);
        title('abs(differences), no reg');
        im(9, real(recon(2).im(:,:,ii)), clim);
        title('real of recon image, Tikhonov');
        im(10, imag(recon(2).im(:,:,ii)), clim);
        title('imag of recon image, Tikhonov');
        im(11, abs(real(image)-real(recon(2).im(:,:,ii)))+ ...
            abs(imag(image)-imag(recon(2).im(:,:,ii))), [-9 9]);
        title('abs(differences), Tikhonov');
    else
        im(1, real(image), clim);
        title('real of original image');
        im(2, imag(image), clim);
        titlef('imag of original image, noise variance: %d', f.sigmas(ii).^2);
        im(5, real(recon(1).im(:,:,ii)), clim);
        title('real of recon image, masked');
        im(6, imag(recon(1).im(:,:,ii)), clim);
        title('imag of recon image, masked');
        im(7, abs(real(image)-real(recon(1).im(:,:,ii)))+ ...
            abs(imag(image)-imag(recon(1).im(:,:,ii))), [-9 9]);
        title('abs(differences), masked');
    end

prompt
end

author: JeffFessler

repo: mirt


text:
classdef SecondOrderGravOnlyPropagator < AbstractPropagator
    %ForceModelPropagator Summary of this class goes here
    %   Detailed explanation goes here

    properties
        forceModels ForceModelsEnum = [ForceModelsEnum.Gravity];
    end

    properties(Constant)
        propagatorEnum = PropagatorEnum.SecOrdGravOnly;
    end

    methods
        function obj = SecondOrderGravOnlyPropagator()

        end

        function [t,y,te,ye,ie] = propagate(obj, integrator, tspan, eventInitStateLogEntry, ...
                                            eventTermCondFuncHandle, termCondDir, maxT, checkForSoITrans, nonSeqTermConds, nonSeqTermCauses, minAltitude, celBodyData, ...
                                            tStartPropTime, maxPropTime)

            if(not(isa(integrator, 'AbstractSecondOrderIntegrator')))
                error('The selected integrator must be a second order integrator in order to use this propagator.');
            end

            plugins = eventInitStateLogEntry.lvdData.plugins;

            %Create function handles
            odefun = obj.getOdeFunctionHandle(eventInitStateLogEntry);
            evtsFunc = obj.getOdeEventsFunctionHandle(eventInitStateLogEntry, eventTermCondFuncHandle, termCondDir, maxT, checkForSoITrans, nonSeqTermConds, nonSeqTermCauses, minAltitude, celBodyData);
            odeOutputFun = obj.getOdeOutputFunctionHandle(tStartPropTime, maxPropTime, eventInitStateLogEntry, plugins);

            %Propagate!
            [t0,y0,yp0] = eventInitStateLogEntry.getSecondOrderIntegratorStateRepresentation();

            if(eventInitStateLogEntry.isHoldDownEnabled())
                %Integrate in the body-fixed frame with zero rates
                %For performance reasons
                bodyInfo = eventInitStateLogEntry.centralBody;
                [rVectECEF, vVectECEF] = getFixedFrameVectFromInertialVect(t0, y0(1:3)', bodyInfo, yp0(1:3)');
                y0 = [rVectECEF', vVectECEF'];

                [t,y,yp,te,ye,ype,ie] = integrator.integrate(odefun, tspan, y0, yp0, evtsFunc, odeOutputFun);

                [rVectECI, vVectECI] = getInertialVectFromFixedFrameVect(t, y(:,1:3)', bodyInfo, yp(:,1:3)');
                y = [rVectECI', vVectECI'];

                if(~isempty(ye))
                    [rVectECIe, vVectECIe] = getInertialVectFromFixedFrameVect(te, ye(:,1:3)', bodyInfo, ype(:,1:3)');
                    ye = [rVectECIe', vVectECIe'];
                end
            else
                [t, y, yp, te, ye, ype, ie] = integrator.integrate(odefun, tspan, y0, yp0, evtsFunc, odeOutputFun);
            end

            y = horzcat(y,yp);
            ye = horzcat(ye,ype);
        end

        function odeFH = getOdeFunctionHandle(obj, eventInitStateLogEntry)
            tankStates = eventInitStateLogEntry.getAllActiveTankStates();
            dryMass = eventInitStateLogEntry.getTotalVehicleDryMass();
            pwrStorageStates = eventInitStateLogEntry.getAllActivePwrStorageStates();
            odeFH = @(t,y) SecondOrderGravOnlyPropagator.odefun(t,y, eventInitStateLogEntry, tankStates, dryMass, pwrStorageStates, obj.forceModels);
        end

        function odeEventsFH = getOdeEventsFunctionHandle(~, eventInitStateLogEntry, eventTermCondFuncHandle, termCondDir, maxT, checkForSoITrans, nonSeqTermConds, nonSeqTermCauses, minAltitude, celBodyData)
            odeEventsFH = @(t,y,yp) AbstractPropagator.odeEvents(t,vertcat(y,yp), eventInitStateLogEntry, eventTermCondFuncHandle, termCondDir, maxT, checkForSoITrans, nonSeqTermConds, nonSeqTermCauses, minAltitude, celBodyData);
        end

        function odeOutputFH = getOdeOutputFunctionHandle(~, tStartPropTime, maxPropTime, eventInitStateLogEntry, plugins)
            odeOutputFH = @(t,y,yp,flag) SecondOrderGravOnlyPropagator.odeOutput(t,y,yp,flag, tStartPropTime, maxPropTime, eventInitStateLogEntry, plugins);
        end

        function [value,isterminal,direction,causes] = callEventsFcn(obj, odeEventsFun, stateLogEntry)
            [t,y,yp] = stateLogEntry.getSecondOrderIntegratorStateRepresentation();
            [value,isterminal,direction,causes] = odeEventsFun(t,y,yp);
        end

        function openOptionsDialog(obj)
            fms = obj.forceModels;

            fmArr = ForceModelsEnum.getEnumsOfDisablableForceModels();
            fmArr = fmArr([fmArr.allowedForSecondOrder] == true);

            [~,initSelInds] = ismember(fms, fmArr);
            initSelInds = initSelInds(initSelInds > 0);

            out = AppDesignerGUIOutput();
            listdlgARH_App('ListString',{fmArr.name}, ...
                            'SelectionMode', 'multiple', ...
                            'ListSize', [300, 300], ...
                            'Name', 'Select Force Models', ...
                            'PromptString', {'Select the Force Models you wish to have enabled during this','event.  Gravity is always enabled.  Disabling Thrust during','periods of coasting may improve performance considerably.'}, ...
                            'InitialValue', initSelInds, ...
                            'out',out);
            Selection = out.output{1};
            ok = out.output{2};

            if(ok == 1)
                obj.forceModels = [ForceModelsEnum.getAllForceModelsThatCannotBeDisabled(), fmArr(Selection)'];
            end
        end

        function tf = canProduceThrust(obj)
            tf = false;
        end
    end

    methods(Static)
        function [ut, rVect] = decomposeIntegratorTandY(t,y)
            ut = t;
            rVect = y(1:3);
        end
    end

    methods(Static, Access=private)
        %%%
        %ODE Function
        %%%
        function d2ydt2 = odefun(t,y, eventInitStateLogEntry, tankStates, dryMass, powerStorageStates, fmEnums)
            bodyInfo = eventInitStateLogEntry.centralBody;
            if(isstruct(bodyInfo.celBodyData) || isempty(bodyInfo.celBodyData))
                bodyInfo.celBodyData = eventInitStateLogEntry.celBodyData;
            end

            [ut, rVect] = SecondOrderGravOnlyPropagator.decomposeIntegratorTandY(t,y);
            vVect = [0;0;0]; %placeholder - this ODE function can't be a function of velocity, only position
            altitude = norm(rVect) - bodyInfo.radius;

            holdDownEnabled = eventInitStateLogEntry.isHoldDownEnabled();

            d2ydt2 = zeros(length(y),1);
            if(holdDownEnabled)
                %launch clamp is enabled, only motion is circular motion
                %(fixed to body)
                %In this case, we are integrating in the body-fixed frame,
                %so all rates are effectively zero
                d2ydt2(1:3) = [0;0;0];
            else
                %launch clamp disabled, propagate like normal
                if(altitude <= 0 && any(fmEnums == ForceModelsEnum.Normal))
                    rswVVect = rotateVectorFromEciToRsw(vVect, rVect, vVect);
                    rswVVect(1) = 0; %kill vertical velocity because we don't want to go throught the surface of the planet
                    vVect = rotateVectorFromRsw2Eci(rswVVect, rVect, vVect);
                end

                thirdBodyGravity = eventInitStateLogEntry.thirdBodyGravity;

                totalMass = eventInitStateLogEntry.getTotalVehicleMass(); %this isn't the total mass but because we can't

                if(totalMass > 0)
                    [forceSum] = TotalForceModel.getForce(fmEnums, ut, rVect, vVect, totalMass, bodyInfo, [], [], [], [], [], [], dryMass, [], thirdBodyGravity, [], []);
                    accelVect = forceSum/totalMass;
                else
                    accelVect = zeros(3,1);
                end

                d2ydt2(1:3) = accelVect;
            end
        end

        %%%
        %ODE Output
        %%%
        function status = odeOutput(t,y,yp,flag, intStartTime, maxIntegrationDuration, eventInitStateLogEntry, plugins)
            y = vertcat(y,yp);
            plugins.executePluginsAfterTimeStepOdeOutputFcn(t,y,flag, eventInitStateLogEntry);

            integrationDuration = toc(intStartTime);

            status = 0;
            if(integrationDuration > maxIntegrationDuration)
                status = 1;
            end
        end
    end
end


text:
function [s] = fslHwrite(P,DIM,VOX,SCALE,TYPE,OFFSET,ORIGIN,DESCRIP)
% writes a header
% FORMAT [s] = spm_hwrite(P,DIM,VOX,SCALE,TYPE,OFFSET,ORIGIN,DESCRIP);
%
% P       - filename 	     (e.g 'spm' or 'spm.img')
% DIM     - image size       [i j k [l]] (voxels)
% VOX     - voxel size       [x y z [t]] (mm [sec])
% SCALE   - scale factor
% TYPE    - datatype (integer - see spm_type)
% OFFSET  - offset (bytes)
% ORIGIN  - [i j k] of origin  (default = [0 0 0])
% DESCRIP - description string (default = 'spm compatible')
%
% s       - number of elements successfully written (should be 348)
%___________________________________________________________________________
%
% spm_hwrite writes variables from working memory into a SPM/ANALYZE
% compatible header file.  The 'originator' field of the ANALYZE format has
% been changed to ORIGIN in the SPM version of the header. funused1
% of the ANALYZE format is used for SCALE
%
% see also dbh.h (ANALYZE) spm_hread.m and spm_type.m
%
%__________________________________________________________________________
% @(#)spm_hwrite.m	2.2 99/10/29

% ensure correct suffix {.hdr} and open header file
%---------------------------------------------------------------------------
P               = P(P ~= ' ');
q    		= length(P);
if q>=4 & P(q - 3) == '.', P = P(1:(q - 4)); end;
P     		= [P '.hdr'];

% For byte swapped data-types, also swap the bytes around in the headers.
mach = 'native';
%if spm_type(TYPE,'swapped'),
%	if spm_platform('bigend'),
		mach = 'ieee-le';
%	else,
%		mach = 'ieee-be';
%	end;
%	TYPE = spm_type(spm_type(TYPE));
%end;
fid             = fopen(P,'w',mach);

if (fid == -1),
	error(['Error opening ' P '. Check that you have write permission.']);
end;
%---------------------------------------------------------------------------
data_type 	= ['dsr      ' 0];

P     		= [P '                  '];
db_name		= [P(1:17) 0];

% set header variables
%---------------------------------------------------------------------------
DIM		= DIM(:)'; if size(DIM,2) < 4; DIM = [DIM 1]; end
VOX		= VOX(:)'; if size(VOX,2) < 4; VOX = [VOX 0]; end
dim		= [4 DIM(1:4) 0 0 0];
pixdim		= [0 VOX(1:4) 0 0 0];
vox_offset      = OFFSET;
funused1	= SCALE;
glmax		= 1;
glmin		= 0;
bitpix 		= 0;
descrip         = zeros(1,80);
aux_file        = ['none                   ' 0];
origin          = [0 0 0 0 0];

%---------------------------------------------------------------------------
if TYPE == 1;   bitpix = 1;  glmax = 1;        glmin = 0;	end
if TYPE == 2;   bitpix = 8;  glmax = 255;      glmin = 0;	end
if TYPE == 4;   bitpix = 16; glmax = 32767;    glmin = 0;  	end
if TYPE == 8;   bitpix = 32; glmax = (2^31-1); glmin = 0;	end
if TYPE == 16;  bitpix = 32; glmax = 1;        glmin = 0;	end
if TYPE == 64;  bitpix = 64; glmax = 1;        glmin = 0;	end

%---------------------------------------------------------------------------
if nargin >= 7; origin = [ORIGIN(:)' 0 0];  end
if nargin <  8; DESCRIP = 'spm compatible'; end

d          	= 1:min([length(DESCRIP) 79]);
descrip(d) 	= DESCRIP(d);

fseek(fid,0,'bof');

% write (struct) header_key
%---------------------------------------------------------------------------
fwrite(fid,348,		'int32');
fwrite(fid,data_type,	'char' );
fwrite(fid,db_name,	'char' );
fwrite(fid,0,		'int32');
fwrite(fid,0,		'int16');
fwrite(fid,'r',		'char' );
fwrite(fid,'0',		'char' );

% write (struct) image_dimension
%---------------------------------------------------------------------------
fseek(fid,40,'bof');

fwrite(fid,dim,		'int16');
fwrite(fid,'mm',	'char' );
fwrite(fid,0,		'char' );
fwrite(fid,0,		'char' );

fwrite(fid,zeros(1,8),	'char' );
fwrite(fid,0,		'int16');
fwrite(fid,TYPE,	'int16');
fwrite(fid,bitpix,	'int16');
fwrite(fid,0,		'int16');
fwrite(fid,pixdim,	'float');
fwrite(fid,vox_offset,	'float');
fwrite(fid,funused1,	'float');
fwrite(fid,0,		'float');
fwrite(fid,0,		'float');
fwrite(fid,0,		'float');
fwrite(fid,0,		'float');
fwrite(fid,0,		'int32');
fwrite(fid,0,		'int32');
fwrite(fid,glmax,	'int32');
fwrite(fid,glmin,	'int32');

% write (struct) image_dimension
%---------------------------------------------------------------------------
fwrite(fid,descrip,	'char');
fwrite(fid,aux_file,    'char');
fwrite(fid,0,           'char');
fwrite(fid,origin,      'int16');
if fwrite(fid,zeros(1,85), 'char')~=85
	fclose(fid);
%	spm_unlink(P);
	error(['Error writing ' P '. Check your disk space.']);
end

s   = ftell(fid);
fclose(fid);


text:
%kkbitxor 'Output = Input 1 XOR (Input 2 or Constant)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kbitxor.pane file
%
% Parameters:
% InputFile: i1 'Input 1', required: 'First input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
% InputFile: i2 'Input 2', optional: 'Second input data object'
% Integer: real 'Constant', default: 0: 'Constant value'
%
% Example: o = kkbitxor({i1, i2}, {'i1','';'o','';'i2','';'real',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kbitxor - Output = Input 1  XOR  (Input 2 or Constant)
%
%  DESCRIPTION
% The bitwise XOR operator performs a bitwise Exclusive OR between each point in
% \fBInput 1" and either \fBInput 2\fP or the \fBConstant\fP value, which
% ever is specified by the user.
%
%  "Data Type - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "Data Type - Two Input Objects" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_2input
% The bitwise operators will not operate on float or complex data types.
%
%  "Map Data" 5
% The bitwise operations have not been written to be fully polymorphic yet.
% They does not
% check for map data, and will therefore always operate on the value data,
% even if a map exists.  This will most likely corrupt indexing into the map.
% In the case of a single input, it is recommended to use the "Copy to
% Value" (kcptoval) segment operator to temporarily move the map data into
% the value segment, run the bitwise operation on the data, then move it
% back to the map with the "Copy from Value" (kcpfromval) operator.
% When operating with two input objects, where at least one of the objects
% contains map data, data should be mapped prior to running the bitwise
% operator.  The "Map Data" operator (kmapdata) can be used to perform
% the mapping operation.
%
%  "Validity Mask - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "Validity Mask - Two Input Objects" 5
% Masking has not been implemented yet for the bitwise operators.  Therefore,
% only the mask from the first input object will be transferred to the output.
%
%  "Input Objects of Different Sizes" 5
% .cI $DATAMANIP/repos/shared/man/sections/resize_2input
% The value used to pad the data when the input files are not the same size is
% zero.
%
%  "Explicit Location and Time Data - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "Explicit Location and Time Data - Two Input Objects" 5
% The bitwise operations have not been extended to understand location and
% time data.  Therefore, only location and time data present in the first input
% object will be transferred to the output.
%
% Executing \fBXOR" runs the program "kbitwise\fP with the -xor flag.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::kbitwise
%
%  RESTRICTIONS
%
% THIS ROUTINE is still under development, and has not been modified to fully
% support the polymorphic data model.  See paragraphs above for discussions
% concerning the polymorphic data segments.
%
% All processing is currently performed as UNSIGNED long, so the operations
% may not work properly on negative values.
%
% The bitwise operators will not operate on float or complex data types.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%% Etract IDE features
clear;clc;
addpath('../caffe/matlab/');
addpath(genpath('utils/'));
% load model and creat network
caffe.set_device(0);
caffe.set_mode_gpu();
netname = 'ResNet_50'; % network: CaffeNet  or ResNet_50

% set your path to the prototxt and model
model =  ['../models/duke/' netname '/' netname '_test.prototxt'];
weights = ['../output/duke_train/IDE_' netname '.caffemodel'];
net = caffe.Net(model, weights, 'test');

if strcmp(netname, 'CaffeNet')
    im_size = 227;
    feat_dim = 4096;
else
    im_size = 224;
    feat_dim = 2048;
end

% mean data
mean_data = importdata('../caffe/matlab/+caffe/imagenet/ilsvrc_2012_mean.mat');
image_mean = mean_data;
off = floor((size(image_mean,1) - im_size)/2)+1;
image_mean = image_mean(off:off+im_size-1, off:off+im_size-1, :);

ef_path = {'dataset/bounding_box_test/', 'dataset/query/'};
ef_name = {'test', 'query'};

if ~exist('feat')
    mkdir('feat')
end

% extract features
for i = 1:2
    img_path = ef_path{i};
    img_file = dir([img_path '*.jpg']);
    feat = single(zeros(feat_dim, length(img_file)));

    for n = 1:length(img_file)
        if mod(n, 1000) ==0
            fprintf('%s: %d/%d\n',ef_name{i}, n, length(img_file))
        end
        img_name = [img_path  img_file(n).name];
        im = imread(img_name);
        im = prepare_img( im, image_mean, im_size);
        feat_img = net.forward({im});
        feat(:, n) = single(feat_img{1}(:));
    end

    save(['feat/IDE_'  netname  '_' ef_name{i} '.mat'], 'feat');
    feat = [];
end

caffe.reset_all();

repo: Learning-via-Translation


text:
function dn = dudiv(dn0,dn1)

% DDIV dual number division
%
%   DN = DDIV(DN0,DN1) returns the dual number DN which is the dual number
%     division DN0/DN1 of the dual numbers DN0 and DN1
%      - DN0 (resp. DN1) is a dual number (DN0 = a0 + eps*b0, eps^2 = 0).
%         It is a 2-vector or a 2*N array (column i represents dual number
%         i) where N is the number of dual numbers. DN0 and DN1 must have
%         the same size. The non-dual part of DN1 must be different from 0.
%      - DN is a dual number. It is a 2*N array (each column is the dual
%          multiplication of the corresponding columns in DN0 and DN1)

s0 = size(dn0);
s1 = size(dn1);
if s0 == [1 2]
    dn0 = dn0.';
    s0 = size(dn0);
end
if s1 == [1 2]
    dn1 = dn1.';
    s1 = size(dn1);
end

% wrong format
if s0(1) ~= 2 || s1(1) ~= 2
    error('DualQuaternion:Ddiv:wrongsize',...
        '%d rows in the DN0 array and %d rows in the DN1 array. It should be 2 for both.',...
        s0(1),s1(1));
end

% sizes do not match
n1 = s0(2);
n2 = s1(2);
if n1 ~= n2
    error('DualQuaternion:Ddiv:notMatch',...
        '%d dual numbers in DN0 array and %d dual numbers in DN1 array.They should be equal.',...
        n1,n2);
end

dn = sym(zeros(2,n1));
dn(1,:) = dn0(1,:)./dn1(1,:);
dn(2,:) = dn0(2,:)./dn1(1,:)-dn0(1,:).*dn1(2,:)./(dn1(1,:).^2);

author: Sable

repo: mcbench-benchmarks


text:
function y = efilter2(x, f, extmod, shift)
% EFILTER2   2D Filtering with edge handling (via extension)
%
%	y = efilter2(x, f, [extmod], [shift])
%
% Input:
%	x:	input image
%	f:	2D filter
%	extmod:	[optional] extension mode (default is 'per')
%	shift:	[optional] specify the window over which the
%		convolution occurs. By default shift = [0; 0].
%
% Output:
%	y:	filtered image that has:
%		Y(z1,z2) = X(z1,z2)*F(z1,z2)*z1^shift(1)*z2^shift(2)
%
% Note:
%	The origin of filter f is assumed to be floor(size(f)/2) + 1.
%	Amount of shift should be no more than floor((size(f)-1)/2).
%	The output image has the same size with the input image.
%
% See also:	EXTEND2, SEFILTER2

if ~exist('extmod', 'var')
    extmod = 'per';
end

if ~exist('shift', 'var')
    shift = [0; 0];
end

% Periodized extension
sf = (size(f) - 1) / 2;

xext = extend2(x, floor(sf(1)) + shift(1), ceil(sf(1)) - shift(1), ...
	       floor(sf(2)) + shift(2), ceil(sf(2)) - shift(2), extmod);

% Convolution and keep the central part that has the size as the input
y = conv2(xext, f, 'valid');

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB