text:
function [density, sound_speed] = hounsfield2density(ct_data, plot_fitting)
%HOUNSFIELD2DENSITY   Convert Hounsfield units to density.
%
% DESCRIPTION:
%       hounsfield2density converts Hounsfield units to units of density
%       [kg/m^3] based on the experimental data given by Schneider et al.
%       The conversion is made using a piece-wise linear fit to the data.
%       For soft-tissue, the approximate sound speed can also be returned
%       using the empirical relationship given by Mast.
%
% USAGE:
%       density = hounsfield2density(ct_data)
%       density = hounsfield2density(ct_data, plot_fitting)
%       [density, sound_speed] = hounsfield2density(ct_data)
%       [density, sound_speed] = hounsfield2density(ct_data, plot_fitting)
%
% INPUTS:
%       ct_data      - CT data in Hounsfield units to convert to density
%
% OPTIONAL INPUTS:
%       plot_fitting - Boolean controlling whether the original data points
%                      and fitting is plotted (default = false)
%
% OUTPUTS:
%       density      - density in kg/m^3
%       sound_speed  - sound speed in m/s
%
% ABOUT:
%       author       - Bradley Treeby
%       date         - 9th January 2012
%       last update  - 10th January 2012
%
% REFERENCES:
%       Schneider, U., Pedroni, E., and Lomax A., "The calibration of CT
%       Hounsfield units for radiotherapy treatment planning," Phys. Med.
%       Biol., 41, pp. 111-124 (1996).
%
%       Mast, T. D., "Empirical relationships between acoustic parameters
%       in human soft tissues," Acoust. Res. Lett. Online, 1(2), pp. 37-42
%       (2000).
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave 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.
%
% k-Wave 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 Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% create empty density matrix
density = zeros(size(ct_data));

% apply conversion in several parts using linear fits to the data
% Part 1: Less than 930 Hounsfield Units
density(ct_data < 930) = polyval([1.025793065681423  -5.680404011488714], ct_data(ct_data < 930));

% Part 2: Between 930 and 1098 (soft tissue region)
density(ct_data >= 930 & ct_data <= 1098) = polyval([0.9082709691264   103.6151457847139], ct_data(ct_data >= 930 & ct_data <= 1098));

% Part 3: Between 1098 and 1260 (between soft tissue and bone)
density(ct_data > 1098 & ct_data < 1260) = polyval([0.5108369316599   539.9977189228704], ct_data(ct_data > 1098 & ct_data < 1260));

% Part 4: Greater than 1260 (bone region)
density(ct_data >= 1260) = polyval([0.6625370912451   348.8555178455294], ct_data(ct_data >= 1260));

% calculate corresponding sound speed values if required using soft tissue
% relationship
if nargout == 2
    sound_speed = (density + 349)./0.893;
end

% plot original data and fitted curves if required
if nargin == 2 && plot_fitting

    % soft tissue values excluding spongiosa
    density_soft_tissue = [0.95, 1.06, 1.04, 1.02, 1.00, 1.07, 1.03, 1.06, 1.05, 1.06, 1.05, 1.03, 1.05, 1.05, 1.04, 1.10, 1.03, 0.98, 1.09, 1.06, 1.04, 1.05]*1000;
    hounsfd_soft_tissue = [ 930, 1055, 1037, 1003, 1003, 1050, 1023, 1055, 1043, 1053, 1044, 1028, 1042, 1045, 1032, 1098, 1014,  958, 1075, 1054, 1032, 1040];

    % bone values
    density_bone = [1.92, 1.61, 1.33, 1.46, 1.68, 1.41, 1.52, 1.29, 1.18, 1.42, 1.33]*1000;
    hounsfd_bone = [2376, 1903, 1499, 1683, 2006, 1595, 1763, 1413, 1260, 1609, 1477];

    % filled lung values
    density_lung = 0.26*1000;
    hounsfd_lung = 259;

    % find linear fit for soft tissue data points
    h_axis_soft_tissue = min(hounsfd_soft_tissue(:)):max(hounsfd_soft_tissue(:));
    p_soft_tissue = polyfit(hounsfd_soft_tissue, density_soft_tissue, 1);

    % find linear fit for bone data points
    h_axis_bone = min(hounsfd_bone(:)):max(hounsfd_bone(:));
    p_bone = polyfit(hounsfd_bone, density_bone, 1);

    % find linear fit from soft tissue to bone region
    h_axis_tissue_to_bone = [h_axis_soft_tissue(end), h_axis_bone(1)];
    density_tissue_to_bone = [polyval(p_soft_tissue, h_axis_soft_tissue(end)), polyval(p_bone, h_axis_bone(1))];
    p_tissue_to_bone = polyfit(h_axis_tissue_to_bone, density_tissue_to_bone, 1);

    % find linear fit from filled lung to soft tissue region
    h_axis_lung_to_tissue = [hounsfd_lung, h_axis_soft_tissue(1)];
    density_lung_to_tissue = [density_lung, polyval(p_soft_tissue, h_axis_soft_tissue(1))];
    p_lung_to_tissue = polyfit(h_axis_lung_to_tissue, density_lung_to_tissue, 1);

    % plot original data points
    figure;
    hold on;
    plot(hounsfd_soft_tissue, density_soft_tissue, 'r.');
    plot(hounsfd_bone , density_bone , 'b.');
    plot(hounsfd_lung, density_lung, 'g.');

    % plot fitting data
    plot(h_axis_soft_tissue, polyval(p_soft_tissue, h_axis_soft_tissue), 'r-');
    plot(h_axis_bone, polyval(p_bone, h_axis_bone), 'b-');
    plot(h_axis_tissue_to_bone, polyval(p_tissue_to_bone, h_axis_tissue_to_bone), 'k-');
    plot(h_axis_lung_to_tissue, polyval(p_lung_to_tissue, h_axis_lung_to_tissue), 'k-');
    xlabel('Hounsfield Units');
    ylabel('Density [kg/m^3]');
    legend('Soft Tissue', 'Bone', 'Lung', 'Location', 'NorthWest');
    box on;

end

repo: Aero-matlab


text:
% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% Our SLAM system requires a number of components:
% * a vehicle
% * a map that defines the positions of some known landmarks in the world
% * a sensor, a range-bearing sensor in this case
% * a SLAM filter

% Creating the vehicle.  First we define the covariance of the vehicles's odometry
% which reports distance travelled and change in heading angle

V = diag([0.005, 0.5*pi/180].^2);

% then use this to create an instance of a Vehicle class
veh = Vehicle(V);

% and then add a "driver" to move it between random waypoints in a square
% region with dimensions from -10 to +10

veh.add_driver( RandomPath(10) );

% Creating the map.  The map covers a square region with dimensions from
% -10 to +10 and contains 20 randomly placed landmarks
map = Map(20, 10);

% Creating the sensor.  We firstly define the covariance of the sensor measurements
% which report distance and bearing angle
W = diag([0.1, 1*pi/180].^2);

% and then use this to create an instance of the Sensor class.
sensor = RangeBearingSensor(veh, map, W, 'animate');
% Note that the sensor is mounted on the moving robot and observes the features
% in the world so it is connected to the already created Vehicle and Map objects.

% Create the filter.  First we need to determine the initial covariance of the
% vehicle, this is our uncertainty about its pose (x, y, theta)
P0 = diag([0.005, 0.005, 0.001].^2);

% Now we create an instance of the EKF filter class
ekf = EKF(veh, V, P0, sensor, W, []);
% and connect it to the vehicle and the sensor and give estimates of the vehicle
% and sensor covariance (we never know this is practice).

% Now we will run the filter for 1000 time steps.  At each step the vehicle
% moves, reports its odometry and the sensor measurements and the filter updates
% its estimate of the vehicle's pose
ekf.run(1000);
% all the results of the simulation are stored within the EKF object

% First let's plot the map
clf; map.plot()
% and then overlay the path actually taken by the vehicle
veh.plot_xy('b');
% and then overlay the path estimated by the filter
ekf.plot_xy('r');
% which we see are pretty close

% Now let's plot the error in estimating the pose
ekf.plot_error()
% and this is overlaid with the estimated covariance of the error.

% Remember that the SLAM filter has not only estimated the robot's pose, it has
% simultaneously estimated the positions of the landmarks as well.  How well did it
% do at that task?  We will show the landmarks in the map again
map.plot();
% and this time overlay the estimated landmark (with a +) and the 3sigma
% uncertainty bounds as green ellipses
ekf.plot_map(3,'g');

author: Allopart

repo: rbpf-gmapping


text:
function selected = spm_eeg_select_event_ui(event)
% Allow the user to select an event using GUI
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_select_event_ui.m 5592 2013-07-24 16:25:55Z vladimir $

selected={};

if isempty(event)
    fprintf('no events were found\n');
    return
end

eventtype = unique({event.type});
Neventtype = length(eventtype);

% Two lists are built in parallel
settings={}; % The list of actual values to be used later
strsettings={}; % The list of strings to show in the GUI

for i=1:Neventtype
    sel = find(strcmp(eventtype{i}, {event.type}));

    numind = find(...
        cellfun('isclass', {event(sel).value}, 'double') & ...
        ~cellfun('isempty', {event(sel).value}));

    charind = find(cellfun('isclass', {event(sel).value}, 'char'));

    emptyind = find(cellfun('isempty', {event(sel).value}));

    if ~isempty(numind)
        numvalue = unique([event(sel(numind)).value]);
        for j=1:length(numvalue)
            ninstances = sum([event(sel(numind)).value] == numvalue(j));
            strsettings=[strsettings; {['Type: ' eventtype{i} ' ; Value: ' num2str(numvalue(j)) ...
                ' ; ' num2str(ninstances) ' instances']}];
            settings=[settings; [eventtype(i), {numvalue(j)}]];
        end
    end

    if ~isempty(charind)
        charvalue = unique({event(sel(charind)).value});
        if ~iscell(charvalue)
            charvalue = {charvalue};
        end
        for j=1:length(charvalue)
            ninstances = length(strmatch(charvalue{j}, {event(sel(charind)).value}, 'exact'));
            strsettings=[strsettings; {['Type: ' eventtype{i} ' ; Value: ' charvalue{j}...
                ' ; ' num2str(ninstances) ' instances']}];
            settings=[settings; [eventtype(i), charvalue(j)]];
        end
    end

    if ~isempty(emptyind)
        strsettings=[strsettings; {['Type: ' eventtype{i} ' ; Value: ; ' ...
            num2str(length(emptyind)) ' instances']}];
        settings=[settings; [eventtype(i), {[]}]];
    end
end

[selection,ok]= listdlg('ListString',strsettings, 'SelectionMode', 'multiple', 'Name', 'Select event', 'ListSize', [400 300]);

if ok
    selected=settings(selection, :);
else
    selected={};
end

author: spm


text:
function spectrumButtonCallback(blockHandle)
% Callback function for SpecButton
% Allows user to specify what .mat file to use for spectrumImport wave type

    [filename,filepath] = uigetfile('.mat');               % Return filename, filepath of user-chosen file
    % Don't set value if no file is chosen, or prompt canceled.
    if ~isequal(filename,0) && ~isequal(filepath,0)
        mask = Simulink.Mask.get(blockHandle);
        fileParam = mask.getParameter('spectrumFile');
        fileParam.Value = [filepath,filename];             % Update spectrumFile with new filename
    end
end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function [Vlims, dVlims] = opf_vlim_fcn(x, mpc, idx, mpopt)
%OPF_VLIM_FCN  Evaluates voltage magnitudes and their gradients.
%   [Vlims, dVlims] = OPF_VLIM_FCN(X, MPC, IDX, MPOPT)
%
%   Computes the voltage magnitudes using real and imaginary part of complex voltage for
%   AC optimal power flow. Computes constraint vectors and their gradients.
%
%   Inputs:
%     X : optimization vector
%     MPC : MATPOWER case struct
%     IDX : index of buses whose voltage magnitudes should be fixed
%     MPOPT : MATPOWER options struct
%
%   Outputs:
%     VLIMS  : vector of voltage magnitudes
%     DVLIMS : (optional) magnitude gradients
%
%   Examples:
%       Vlims = opf_vlim_fcn(x, mpc, mpopt);
%       [Vlims, dVlims] = opf_vlim_fcn(x, mpc, idx, mpopt);
%
%   See also OPF_VLIM_HESS

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

%% define named indices into data matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;

%% unpack data
[Vr, Vi] = deal(x{:});

%% problem dimensions
nb = length(Vi);            %% number of buses
n = length(idx);            %% number of buses with voltage limits

%% compute voltage magnitude
Vm2 = Vr(idx).^2 + Vi(idx).^2;
Vlims = [ mpc.bus(idx, VMIN).^2 - Vm2;
          Vm2 - mpc.bus(idx, VMAX).^2 ];

if nargout > 1
    %% compute partials of voltage magnitude w.r.t Vr and Vi
    dVm_dVr = sparse(1:n, idx, 2 * Vr(idx), n, nb);
    dVm_dVi = sparse(1:n, idx, 2 * Vi(idx), n, nb);
    dVlims = [ -dVm_dVr -dVm_dVi;   %% Vlims w.r.t Vr, Vi
                dVm_dVr  dVm_dVi  ];
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [dfuncGrad] = gpsimMapFunctionalLikeGrad2(model, p)

% GPSIMMAPFUNCTIONALLIKEGRAD2 Compute the functional gradient for GPSIMMAP.
% FORMAT
% DESC computes the functional gradient of the log likelihood for
% use in the MAP approximation to the GPSIM posterior solution.
% ARG model : the model for which the gradient is to be computed.
% RETURN p: the gradient of the log likelihood with respect to the
% points of the function.
%
% SEEALSO : gpsimMapLikeGradientImplicit
%
% COPYRIGHT : Pei Gao, Magnus Rattray and Neil D. Lawrence, 2008

% SHEFFIELDML

dfuncGrad = zeros(1, model.numParams-model.kern.nParams);
gB = zeros(model.numGenes, 1);
gD = zeros(model.numGenes, 1);
gS = zeros(model.numGenes, 1);
ngParamk = model.ngParam/model.numGenes;

if isfield(model,'includeNoise') && model.includeNoise
  gn = zeros(model.numGenes, 1);
  noiseMat = ones(length(model.t),1)*model.noiseVar;
  yvar = model.yvar + noiseMat;
else
  yvar = model.yvar;
end

if isfield(model, 'includeRepression') && model.includeRepression
  galpha = zeros(model.numGenes, 1);
  nBaseParams = 4;
else
  nBaseParams = 3;
end

if model.ngParam > 0
  ggParam= zeros(model.numGenes, ngParamk);
end

numData = length(model.t);
for k = 1:model.numGenes
  if model.ngParam
    gInd = k;
  else
    gInd = 1;
  end
  for i=1:numData
    arg = model.t(i)-model.mapt(p);
    if arg >= 0
      ind = i + (k-1)*numData;
      beta_ik=1/yvar(ind);

      [dxdB dxdD dxdS dxdalpha dxdgParam] = gpsimXGradient(model, i, k);

      gB(k) = gB(k) - beta_ik*model.g_grad(p,gInd)*dxdB*exp(-model.D(k)*arg ...
                               +log(model.S(k)) +log(model.step));

      factor = model.ypred(model.times_index(i), k) - model.y(ind);
      gD(k) = gD(k) - beta_ik*model.g_grad(p,gInd)* (dxdD-arg*factor) ...
              *exp(-model.D(k)*arg+log(model.S(k)) +log(model.step));

      gS(k) = gS(k) - beta_ik*model.g_grad(p,gInd)* (dxdS*model.S(k) + ...
                       factor) *exp(-model.D(k)*arg+log(model.step));

      if isfield(model, 'includeRepression') && model.includeRepression
        galpha(k) = galpha(k) - beta_ik*model.g_grad(p,gInd)*dxdalpha ...
            *exp(-model.D(k)*arg+log(model.S(k)) +log(model.step));
      end

      if model.ngParam > 0
        ggParam(k,:) = ggParam(k,:) - beta_ik* (model.g_grad(p,gInd)* ...
            dxdgParam+factor*model.dggrad(p, gInd)) *exp(-model.D(k)* ...
                                arg+log(model.S(k))+log(model.step));
      end

      if isfield(model,'includeNoise') && model.includeNoise
        gn(k) = gn(k) + 2*sqrt(model.noiseVar(k))*beta_ik^2* ...
                model.g_grad(p,gInd)*model.S(k)*factor*exp(-model.D(k)*arg+log(model.step));
      end
    end
  end
end

geneParam = nBaseParams + ngParamk;
startPoint = 1;
endPoint = geneParam;

for k = 1:model.numGenes
  dfuncGrad(startPoint:(startPoint+2)) = [gB(k) gS(k) gD(k)];

  if isfield(model, 'includeRepression') && model.includeRepression
    dfuncGrad(startPoint+3) = galpha(k);
  end

  if model.ngParam > 0
    dfuncGrad((startPoint+nBaseParams):endPoint) = ggParam(k,:)';
  end
  startPoint = startPoint + geneParam;
  endPoint = endPoint + geneParam;
end

if isfield(model,'includeNoise') && model.includeNoise
  dfuncGrad(end-length(gn)+1:end) = gn';
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [p,q] = bandwidth(A)
%BANDWIDTH    Upper and lower bandwidth of matrix A
%
%   A   = 0  for  i-j > p  or j-i > q
%    ij
%
%    [p,q] = bandwidth(A)
%

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

  [p,q] = bandwidth(A.x);

repo: OpenMAS


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
%
% Create random divisions of the training set for cross-validation
%
% ------------------------------------------------------------------------

% Create random divisions of the training set
id = '3';
percentage_train = 0.5;
full_set = 'train2012';

% Get full ids
im_ids = database_ids('pascal2012',full_set);

ids_a = sort(randperm(length(im_ids),floor(percentage_train*length(im_ids))));
ids_b = setdiff(1:length(im_ids),ids_a);
im_ids_a = im_ids(ids_a);
im_ids_b = im_ids(ids_b);

% Write to file
file_a = fullfile(root_dir,'datasets', 'pascal2012','gt_sets',[full_set '_' id 'a.txt']);
file_b = fullfile(root_dir,'datasets', 'pascal2012','gt_sets',[full_set '_' id 'b.txt']);

dlmwrite(file_a,im_ids_a,'')
dlmwrite(file_b,im_ids_b,'')

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
%% FUNCTION nchoose
%   How many unique combinations of 1 to N elements are there. This
%   function is not intended to directly called by users.
%   See the Logistic_iMSF function.
%
%% LICENSE
%   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/>.
%
%   Copyright (C) 2011 - 2012 Lei Yuan, Jiayu Zhou, and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 12, 2012.
%
%% RELATED FUNCTIONS
%   init_opts, Construct_iMSF, Logistic_iMSF, MultiSource_LogisticR

function W = nchoose(S)

N = numel(S) ;

% How many unique combinations of 1 to N elements are there
M = (2^N)-1 ;
if N > 18,
    warning('Nchoose:LargeOutput', ...
        'There are %d unique combinations. Please be patient ...',M) ;
end

S = S(:).' ;   % make the set a row vector, for uniform output

W = cell(M,1) ;    % Pre-allocation of output
p2=2.^(N-1:-1:0) ; % This part of the formula can be taken out of the loop

for i=1:M,
    % calculate the (reversed) binary representation of i
    % select the elements of the set based on this representation
    W{i} = S(bitget(i*p2,N) > 0) ;
end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function [s,r] = medit(varargin)
  % MEDIT wrapper to launch medit to visualize a given tetmesh
  %
  % [s,r] = medit(V,T,F)
  % [s,r] = medit(V,T,F,'ParameterName',ParameterValue)
  %
  % Inputs:
  %   V  #V by 3 list of vertex positions
  %   T  #T by 4|5 list of tetrahedron findices (additional column is color
  %     index)
  %   F  #F by 3|4 list of face indices (additional column is color index)
  %   Optional:
  %     'Wait'  followed by boolean whether to wait for completion {true}
  %     'Data'  followed by #V|#F+#T data values {[]}
  %     'Subplot' follwed by [x y i] defining subplot window (see subplot)
  %     'Edges' follwed by E list of edges to write along with mesh (see
  %       writeMESH)
  %     'Palette'  followed by list (up to 5) iso values
  % Outputs:
  %   s,r  result of system call
  %
  % Example:
  %   % tet mesh in (V,T), per-tet data in D
  %   [F,J] = boundary_faces(T);
  %   % Display data at tets and also on boundary facets
  %   medit(V,T,F,'Data',[D(J);D])
  %
  % See also: tetramesh
  %

  V = varargin{1};
  T = varargin{2};
  if nargin<3
    F = [];
  else
    F = varargin{3};
  end
  wait = true;
  D = [];
  E = [];
  subp = [];
  palette = [];

  ii = 4;
  while(ii<=nargin)
    switch varargin{ii}
    case 'Wait'
      ii = ii + 1;
      assert(ii<=nargin);
      wait = varargin{ii};
    case 'Data'
      ii = ii + 1;
      assert(ii<=nargin);
      D = varargin{ii};
    case 'Edges'
      ii = ii + 1;
      assert(ii<=nargin);
      E = varargin{ii};
    case 'Subplot'
      ii = ii + 1;
      assert(ii<=nargin);
      subp = varargin{ii};
    case 'Palette'
      ii = ii + 1;
      assert(ii<=nargin);
      palette = varargin{ii};
    otherwise
      error(['Unsupported parameter: ' varargin{ii}]);
    end
    ii = ii+1;
  end

  assert(numel(palette)<=5);

  if ~exist('wait','var')
    wait = true;
  end

  % Change these paths accordingly
  MEDIT_PATH = path_to_medit;
  TEMP_MESH_FILE  = '/var/tmp/temp.mesh';
  TEMP_MEDIT_FILE = '/var/tmp/temp.medit';
  TEMP_BB_FILE = '/var/tmp/temp.bb';

  % write default medit options
  f = fopen(TEMP_MEDIT_FILE,'w');
  fprintf(f, [...
    'BackgroundColor 1 1 1\n' ...
    'LineColor 0 0 0\n']);
  if size(F,2) > 3 || size(T,2) > 4
    fprintf(f,['NbMaterials\n']);
    %% collect indices
    %I = [F(:,3:end);T(:,4:end)];
    %count = max(I) - min(I);
    %fprintf(f,'%d\n',count);
    fprintf(f,'RenderMode colorshadinglines\n');
  else
    fprintf(f,'RenderMode shading + lines\n');
  end

  if isempty(subp)
    fprintf(f,'WindowSize 1024 800\n');
  else
    ss = get(0,'ScreenSize');
    w = ss(3)/subp(2);
    h = ss(4)/subp(1);
    [i j] = ind2sub(subp(1:2),subp(3));
    j = (j-1)*w;
    i = (i-1)*h;
    fprintf(f,'WindowSize %d %d\n',round(w),round(h));
    fprintf(f,'WindowPosition %d %d\n',round(j),round(i));
  end

  if ~isempty(palette)
    fprintf(f,'Palette');
    fprintf(f,' %g',palette(:));
    fprintf(f,'\n');
  end

  fclose(f);

  % write bb file
  if isempty(D)
    if exist(TEMP_BB_FILE,'file')
      delete(TEMP_BB_FILE);
    end
  else
    switch size(D,1)
    case size(F,1)
      warning('Appending junk for tet data since medit needs it.');
      D = [D;mean(D)*ones(size(T,1),1)];
      type = 1;
    case size(T,1)
      if size(F,1) == 0
        allF = [ ...
          T(:,2) T(:,4) T(:,3); ...
          T(:,1) T(:,3) T(:,4); ...
          T(:,1) T(:,4) T(:,2); ...
          T(:,1) T(:,2) T(:,3); ...
          ];
        [I,C] = on_boundary(T);
        DF = repmat(D,4,1);
        F = allF(C(:),:);
        DF = DF(C(:),:);
        D = [DF;D];
      else
        warning('Appending junk for face data since medit needs it.');
        D = [mean(D)*ones(size(F,1),1);D];
      end
      type = 1;
    case size(V,1)
      type = 2;
    case size(F,1)+size(T,1)
      type = 1;
    otherwise
      error('size(Data,1) %d should be == size(V,1) (%d) or size(T,1) + size(F,1) (%d + %d = %d)',size(D,1),size(V,1),size(T,1),size(F,1),size(T,1)+size(F,1));
    end
    writeBB(TEMP_BB_FILE,D,type);
  end

  % write temporary mesh
  writeMESH(TEMP_MESH_FILE,V,T,F,E);

  command = [MEDIT_PATH ' ' TEMP_MESH_FILE ' ' TEMP_MEDIT_FILE];
  command
  if ~wait
    command = [command ' &'];
  end
  [s,r] = system(command);
  if(s ~= 0)
      error(r);
  end
end

author: alecjacobson

repo: gptoolbox


text:
function sF = calcDensity(v,varargin)
% calculate a density function out of (weighted) unit vectors
%
% Syntax
%
%   sF = calcDensity(v)
%   sF = calcDensity(v,'weights',w)
%   sF = calcDensity(v,'halfwidth',delta)
%   sF = calcDensity(v,'kernel',psi)
%    f = calcDensity(v,S2G)
%
% Input
%  v   - sampling points for density estimation @vector3d
%  S2G - @vector3d
%  w   - weights, default is all one
%  delta - halfwidth of the kernel, default is 10 degree
%  psi - @S2Kernel function, default is S2 de la Vallee Poussin
%
% Output
%  sF  - @S2Fun
%   f  - function values
%
% Options
%  halfwidth - halfwidth of a kernel
%  kernel    - specify a S2Kernel
%  weights   - vector of weights, with same length as v
%

% determine kernel function
hw = get_option(varargin,'halfwidth',10*degree);
psi = get_option(varargin,'kernel',S2DeLaValleePoussinKernel('halfwidth',hw));

% ignore nans
v = subSet(v,~isnan(v));

sF = 4*pi * S2FunHarmonic.quadrature(v,ones(size(v)),varargin{:});

% normalize
if ~check_option(varargin,'noNormalization')
  sF = sqrt(4*pi) * sF ./ sF.fhat(1);
end

% convolute with kernel function
sF = conv(sF,psi);

% if required compute function values
if nargin > 1 && isa(varargin{1},'vector3d')
  sF = sF.eval(varargin{1});
end

end

author: mtex-toolbox

repo: mtex

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


text:
function [beta,f,y] = trainGM2()

[f,y] = sampleDetector(@detector1,1000000,8);
fprintf(2,'Fitting model...\n');
save 'samples_gm_1_2.mat' f y;
beta = logist2(y',f');
save 'beta_gm_1_2.txt' beta -ascii;

[f,y] = sampleDetector(@detector2,1000000,8);
fprintf(2,'Fitting model...\n');
save 'samples_gm_2_4.mat' f y;
beta = logist2(y',f');
save 'beta_gm_2_4.txt' beta -ascii;

[f,y] = sampleDetector(@detector3,1000000,16);
fprintf(2,'Fitting model...\n');
save 'samples_gm_4_8.mat' f y;
beta = logist2(y',f');
save 'beta_gm_4_8.txt' beta -ascii;

function [f] = detector1(im)
f = detector(im,1);

function [f] = detector2(im)
f = detector(im,2);

function [f] = detector3(im)
f = detector(im,4);

function [f] = detector(im,sigma)
[a] = detGM(im,sigma);
[b] = detGM(im,sigma*2);
a = a(:);
b = b(:);
f = [ ones(size(a)) a b ]';

author: Cloud-CV

repo: object-proposals


text:
function y = iterapp(op,afun,atype,afcnstr,x,varargin)
%ITERAPP   Apply matrix operator to vector and error gracefully.
%   ITERAPP(OP,AFUN,ATYPE,AFCNSTR,X) applies matrix operator AFUN to vector
%   X. If ATYPE is 'matrix, then AFUN is a matrix and the OP is applied
%   directly. OP is either 'mtimes' or 'mldivide'.
%   ATYPE and AFCNSTR are used in case of error.
%   ITERAPP(OP,AFUN,ATYPE,AFCNSTR,X,P1,P2,...) allows extra arguments to
%   AFUN(X,P1,P2,...) although this usage is now discouraged in favor of
%   using anonymous functions.
%   AFUN(X,P1,P2,...,PN,TFLAG) should accept a TFLAG as its final input if
%   the calling function is BICG, LSQR or QMR. TFLAG is either 'transp' or
%   'notransp' depending on whether A' OP X or A OP X is required.
%   ITERAPP is designed for use by iterative methods like PCG which
%   require matrix operators AFUN representing matrices A to operate on
%   vectors X and return A*X and may also have operators MFUN representing
%   preconditioning matrices M operate on vectors X and return M\X.
%
%   See also BICG, BICGSTAB, CGS, GMRES, LSQR, MINRES, PCG, QMR, SYMMLQ.

%   Copyright 1984-2004 The MathWorks, Inc.
%   $Revision: 1.7.4.2 $ $Date: 2004/12/06 16:35:55 $

if isequal(atype,'matrix')
    switch lower(op)
        case 'mtimes'
            if (nargin >= 6) & isequal(varargin{end},'transp')
                y = afun' * x;
            else
                y = afun * x;
            end
        case 'mldivide'
            if (nargin >= 6) & isequal(varargin{end},'transp')
                y = afun' \ x;
            else
                y = afun \ x;
            end
        otherwise
            error('MATLAB:iterapp:InvalidOp', 'Invalid operation.')
    end
else
    try
        if (nargin >= 6) & isequal(varargin{end},'notransp')
            % A request for A*x coming from BICG, LSQR and QMR
            try
                % New syntax: we now request afun(x,P1,P2,...,PN,'notransp')
                y = afun(x,varargin{:});
            catch
                % Old syntax: we used to accept afun(x,P1,P2,...,PN)
                y = afun(x,varargin{1:end-1});
            end
        else
            % A request for A*x
            % coming from BICGSTAB, CGS, GMRES, MINRES, PCG or SYMMLQ
            % with the call always afun(P1,P2,...,PN)
            % or a request for A'*x coming from
            % BICG, LSQR and QMR in the afun(x,P1,P2,...,PN,'transp') case
            y = afun(x,varargin{:});
        end
    catch
        error('MATLAB:InvalidInput', ['user supplied %s ==> %s\n' ...
            'failed with the following error:\n\n%s'], ...
            atype,afcnstr,lasterr);
    end

    if ~isvector(y) || (size(y,2) ~= 1)
        error('MATLAB:MustReturnColumn', ['user supplied %s ==> %s\n' ...
            'must return a column vector.'], ...
            atype,afcnstr)
    end
end

author: johannesgerer

repo: jburkardt-m


text:
function variance = weibull_variance ( a, b, c )

%*****************************************************************************80
%
%% WEIBULL_VARIANCE returns the variance of the Weibull PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, the parameters of the PDF.
%    0.0 < B,
%    0.0 < C.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  g1 = gamma ( ( c + 2.0 ) / c );
  g2 = gamma ( ( c + 1.0 ) / c );

  variance = b * b * ( g1 - g2 * g2 );

  return
end


text:
% StackExchange Signal Processing Q75231
% https://dsp.stackexchange.com/questions/75231
% Image Denoising Using Total Variation (TV) Minimization with ADMM
% References:
%   1.
% Remarks:
%   1.  B
% TODO:
% 	1.  C
% Release Notes
% - 1.0.000     18/05/2021
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

CONV_SHAPE_FULL     = 1;
CONV_SHAPE_SAME     = 2;
CONV_SHAPE_VALID    = 3;

%% Simulation Parameters

imageFileName = 'Lena256.png';

paramLambda = 0.025; %<! TV Norm

numIterations   = 250;
noiseStd        = 10 / 255;

%% Generate Data

mI = im2double(imread(imageFileName));

numRows = size(mI, 1);
numCols = size(mI, 2);

mY = mI + (noiseStd * randn(numRows, numCols));
vY = mY(:);

vXInit  = zeros(numRows * numCols, 1);
% Generate the Diff Operator (2D Gradient) by Finite Differences
mD      = CreateGradientOperator(numRows, numCols);

% Objective Function
hObjFun = @(vX) (0.5 * sum( (vX - vY) .^ 2)) + (paramLambda * sum(abs(mD * vX)));

% Analysis
solverIdx       = 0;
cMethodString   = {};

mObjFunValMse   = zeros(numIterations, 1);
mSolMse         = zeros(numIterations, 1);

%% Display Data

figureIdx = figureIdx + 1;

hFigure = figure('Position', [100, 100, 542, 642]);
hAxes   = axes(hFigure, 'Units', 'pixels', 'Position', [010, 326, 256, 256]);
hImgObj = image(repmat(mI, [1, 1, 3]));
set(get(hAxes, 'Title'), 'String', {['Reference Image']}, ...
    'FontSize', fontSizeTitle);
set(hAxes, 'XTick', [], 'XTickLabel', []);
set(hAxes, 'YTick', [], 'YTickLabel', []);

hAxes   = axes(hFigure, 'Units', 'pixels', 'Position', [276, 326, 256, 256]);
hImgObj = image(repmat(mY, [1, 1, 3]));
set(get(hAxes, 'Title'), 'String', {['Input Image (Noisy)']}, ...
    'FontSize', fontSizeTitle);
set(hAxes, 'XTick', [], 'XTickLabel', []);
set(hAxes, 'YTick', [], 'YTickLabel', []);

%% Solution by CVX

solverString = 'CVX';

% cvx_solver('SDPT3'); %<! Default, Keep numRows low
% cvx_solver('SeDuMi');
cvx_solver('Mosek'); %<! Can handle numRows > 500, Very Good!
% cvx_solver('Gurobi');

hRunTime = tic();

cvx_begin('quiet')
    % cvx_precision('best');
    variable vX(numRows * numCols);
    minimize( (0.5 * sum_square(vX - vY)) + (paramLambda * norm(mD * vX, 1)) );
cvx_end

runTime = toc(hRunTime);

DisplayRunSummary(solverString, hObjFun, vX, runTime, cvx_status);

sCvxSol.vXCvx     = vX(:);
sCvxSol.cvxOptVal = hObjFun(vX);

%% Display Result

hAxes   = axes(hFigure, 'Units', 'pixels', 'Position', [010, 010, 256, 256]);
hImgObj = image(repmat(reshape(vX, numRows, numCols), [1, 1, 3]));
set(get(hAxes, 'Title'), 'String', {['Denoised Image - CVX']}, ...
    'FontSize', fontSizeTitle);
set(hAxes, 'XTick', [], 'XTickLabel', []);
set(hAxes, 'YTick', [], 'YTickLabel', []);

%% Solution by ADMM
%{
Solving:

$$ \arg \min_{ x \in \mathbb{R}^{n} } \frac{1}{2} {\left\| x - y \right|}_{2}^{2} + \lambda {\left\| D x \right\|}_{1} $$
%}

solverIdx                   = solverIdx + 1;
cLegendString{solverIdx}    = ['ADMM'];

hRunTime = tic();

[vX, mX] = SolveProxTvAdmm(vXInit, vY, mD, paramLambda, 'numIterations', numIterations);

runTime = toc(hRunTime);

DisplayRunSummary(cLegendString{solverIdx}, hObjFun, vX, runTime);

[mObjFunValMse, mSolMse] = UpdateAnalysisData(mObjFunValMse, mSolMse, mX, hObjFun, sCvxSol, solverIdx);

%% Display Results

hAxes   = axes(hFigure, 'Units', 'pixels', 'Position', [276, 010, 256, 256]);
hImgObj = image(repmat(reshape(vX, numRows, numCols), [1, 1, 3]));
set(get(hAxes, 'Title'), 'String', {['Denoised Image - ADMM']}, ...
    'FontSize', fontSizeTitle);
set(hAxes, 'XTick', [], 'XTickLabel', []);
set(hAxes, 'YTick', [], 'YTickLabel', []);

if(generateFigures == ON)
    % saveas(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    print(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
end

figureIdx = figureIdx + 1;

hFigure = DisplayComparisonSummary(numIterations, mObjFunValMse, mSolMse, cLegendString, figPosLarge, lineWidthNormal, fontSizeTitle, fontSizeAxis);

if(generateFigures == ON)
    % saveas(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    print(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function [seq_info]=load_seq_info(base_path,seq_name)
%This function loads the information of the video corresponding to the
%given path seq_path, which is modified from
    % GENCONFIG
    % Generate a configuration of a sequence
    %
    % INPUT:
    %   dataset - The name of dataset ('otb','vot2013','vot2014','vot2015')
    %   seqName - The name of a sequence in the given dataset
    %
    % OUTPUT:
    %   config - The configuration of the given sequence
    %
    % Hyeonseob Nam, 2015
    %
% Input:
%       base_path - the base path of the sequences
%       seq_name -  the name of the sequence
% Output:
%       seq_info - a two-filed block including the path of each image file
%       (seq_info.img) and the groundtruth (seq_info.gt).

% By Xin Li, 2017-3-8
%%
seq_info.seq_name=seq_name;

%% load images
% handle special sequences
switch(seq_info.seq_name)
    case {'Jogging-1', 'Jogging-2'}
        seq_info.img_dir = fullfile(base_path, 'Jogging', 'img');
    case {'Skating2-1', 'Skating2-2'}
        seq_info.img_dir = fullfile(base_path, 'Skating2', 'img');
    otherwise
        seq_info.img_dir = fullfile(base_path, seq_info.seq_name, 'img');
end

 if(~exist(seq_info.img_dir,'dir'))
        error('%s does not exist!!',seq_info.img_dir);
 end

 % parse img list
  seq_info.img_list = parseImg(seq_info.img_dir);
  switch(seq_info.seq_name)
      case 'David'
          seq_info.img_list = seq_info.img_list(300:end);
      case 'Tiger1'
          seq_info.img_list = seq_info.img_list(6:end);
  end

  % load gt
  switch(seq_info.seq_name)
      case 'Jogging-1'
          gtPath = fullfile(base_path, 'Jogging', 'groundtruth_rect.1.txt');
      case 'Jogging-2'
          gtPath = fullfile(base_path, 'Jogging', 'groundtruth_rect.2.txt');
      case 'Skating2-1'
          gtPath = fullfile(base_path, 'Skating2', 'groundtruth_rect.1.txt');
      case 'Skating2-2'
          gtPath = fullfile(base_path, 'Skating2', 'groundtruth_rect.2.txt');
      case 'Human4'
          gtPath = fullfile(base_path, 'Human4', 'groundtruth_rect.txt');
      otherwise
          gtPath = fullfile(base_path, seq_info.seq_name, 'groundtruth_rect.txt');
  end

  if(~exist(gtPath,'file'))
      error('%s does not exist!!',gtPath);
  end

  gt = importdata(gtPath);
  switch(seq_info.seq_name)
      case 'Tiger1'
          gt = gt(6:end,:);
      case {'Board','Twinnings'}
          gt = gt(1:end-1,:);
  end

  seq_info.gt = gt;

  nFrames = min(length(seq_info.img_list), size(seq_info.gt,1));
  seq_info.img_list = seq_info.img_list(1:nFrames);
  seq_info.gt = seq_info.gt(1:nFrames,:);
  seq_info.nFrames = nFrames;

author: XinLi-zn

repo: TADT

save_path: github-repos/MATLAB/XinLi-zn-TADT


text:
function out = Dstar(R, n, e, fun)
%finds excess distortion for a given rate at a given blocklength for
%Gaussian source with unit variance
%R - rate
%n - block length (scalar)
%e - excess probability
%fun - which function to use for calculation

%
%   Created in 2012 by Victoria Kostina (vkostina@caltech.edu)
%

%starting points - 'persistent' to make optimization faster
persistent x0;
persistent y0;

try

    switch lower(fun)
        case 'shannon'
            %achievability via Shannon's upper bound:
            tol = 1e-4;
            qglb = [2^(-2*R); 1];
            qgub = [1; Inf];
            options = optimset('TolX',tol, 'MaxFunEvals', 500, 'Algorithm', 'active-set');

            if isempty(x0) && isempty (y0)
                x0 = [2^(-2*R) 1];
                y0 = [(2^(-2*R)+1)/2 1]; %q gamma
            end
            out = DShannon();
        case 'normal'
            %normal approximation
            out = DNormal();
        case 'spherecoveringc'
            %converse via sphere covering
            out = DSphereCovering();
        case 'spherecoveringa'
            %achievability via sphere covering
            if isempty(x0) && isempty (y0)
                x0 = [2^(-2*R) .5]; %distortion
             %   y0 = [0 1];        %auxiliary parameter q
            end
            out = DSphereCoveringA();
        otherwise
            disp('Unknown type.')
    end

catch ME
    fprintf('Error: %s\n', ME.message);
    out = NaN;
end

%--------------------------------------------------------------------------
    function out = DShannon()
        %achievability via Shannon's upper bound:
        out = fzero(@(x)Pexcess(x) - e,x0);
        x0 = out;
        function out = Pexcess(x)
            D = x;
            [y0, Pe] = fmincon(@optqg, y0,[],[],[],[],qglb, qgub, [], options);
            out = Pe;

            function out = optqg(y)
                q = y(1);
                g = log2(y(2));
                out = PeubShannon(R, n, D, q, g);
            end
        end

    end

%--------------------------------------------------------------------------
    function out = DNormal()
        out = 2^(-2*R)*(1 + sqrt(2/n)*Qinv(e));
    end

%--------------------------------------------------------------------------
    function out = DSphereCovering()
        %converse via sphere covering

        %find rstar:
        factor = 10;
        rstar = fzero(@(x)chi2cdf(x^2*n, n) - chi2cdf(factor*n, n) + e, 1);
        out = fminbnd(@(x) abs(R - rate(x)),2^(-2*R),1);

        function out = rate(x)
            out = 1/2*log2(rstar^2/x);
        end
    end

%--------------------------------------------------------------------------
    function out = DSphereCoveringA()
        %achievability via sphere covering
        M = 2.^(n*R);
        out = fzero(@(x)Pexcess(x) - e, x0);
    function out = Pexcess(x)
        d = x;
        s0 = 1; %variance of source
        r0 = sqrt(s0^2 - d); %distance from the center of the representation sphere

        out = quad(@density, r0 - sqrt(d), r0 + sqrt(d))...
            + 1 - chi2cdf((r0 + sqrt(d))^2*n,n)...
            + chi2cdf((r0 - sqrt(d))^2*n,n);
        out = real(out);

        function out = density(r)
          cosa = (r0^2 + r.^2 - d)./(2*r0.*r);
          sina = (1 - cosa.^2).^(1/2);
          out = ub(sina).*chi2pdf(n*r.^2, n).*n*2.*r;
        end

        function out = ub(sina)
            if (n < 50)
                A = 1/sqrt(pi).*gamma(n/2+1)./n./gamma((n-1)/2+1);
                out = (1 - A.*(sina).^(n-1)).^M;
            else
                out = exp(-M*sina.^(n-1)./sqrt(2*pi*n));
            end
        end
    end

    end

end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function [pts, confvals] = confmaps2pts(C)
%CONFMAPS2PTS Convert a set of confidence maps into a set of points.
% Usage:
%   [pts, confvals] = confmaps2pts(C)
%
% See also: pts2confmaps

numChannels = size(C,3);
pts = zeros(numChannels,2,'single');
confvals = zeros(numChannels,1,'like',C);
for i = 1:numChannels
    [confvals(i), ind] = max(vert(C(:,:,i)));
    [r,c] = ind2sub(size(C(:,:,i)),ind);
    pts(i,:) = [c r];
end

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
%% dbDivision.m
% This script randomly divides database into 3 parts - training (TRN),
% validating (VAL) and testing (TST).
%
% 16-07-10 Michal Uricar
% 21-03-12 Michal Uricar, LFW annotation in one file

clearvars; close all;

%% Timestamp

fprintf(1,'Started on %s\n\n', datestr(now));

%% Load pruned database

addpath('./Functions');

% load('./MAT/eyefd_on_lfw_pruned.mat');
load('./MAT/db_good.mat');

tmp_as = annotation_struct;

N = annotation_struct.N;

per60 = round(60/100*N);
per20 = round(20/100*N);

tic
indices_perm = randperm(N);
toc

% 60% of images will form training set
training = indices_perm(1:per60); indices = training;
indices_perm(1:per60) = [];
% save training set
%image = tmp_image(training);
annotation_struct.bbox = tmp_as.bbox(training, :);
annotation_struct.names = tmp_as.names(training);
annotation_struct.eye_r = tmp_as.eye_r(training, :);
annotation_struct.eye_l = tmp_as.eye_l(training, :);
annotation_struct.canthus_rr = tmp_as.canthus_rr(training, :);
annotation_struct.canthus_rl = tmp_as.canthus_rl(training, :);
annotation_struct.canthus_lr = tmp_as.canthus_lr(training, :);
annotation_struct.canthus_ll = tmp_as.canthus_ll(training, :);
annotation_struct.mouth = tmp_as.mouth(training, :);
annotation_struct.mouth_corner_r = tmp_as.mouth_corner_r(training, :);
annotation_struct.mouth_corner_l = tmp_as.mouth_corner_l(training, :);
annotation_struct.nose = tmp_as.nose(training, :);
annotation_struct.N = numel(annotation_struct.names);

% save('./MAT/training.mat', 'image');
save('./MAT/TRN.mat', 'annotation_struct', 'indices');

% 20% of images will form validating set
validating = indices_perm(1:per20); indices = validating;
indices_perm(1:per20) = [];
% save validating set
% image = tmp_image(validating);
annotation_struct.bbox = tmp_as.bbox(validating, :);
annotation_struct.names = tmp_as.names(validating);
annotation_struct.eye_r = tmp_as.eye_r(validating, :);
annotation_struct.eye_l = tmp_as.eye_l(validating, :);
annotation_struct.canthus_rr = tmp_as.canthus_rr(validating, :);
annotation_struct.canthus_rl = tmp_as.canthus_rl(validating, :);
annotation_struct.canthus_lr = tmp_as.canthus_lr(validating, :);
annotation_struct.canthus_ll = tmp_as.canthus_ll(validating, :);
annotation_struct.mouth = tmp_as.mouth(validating, :);
annotation_struct.mouth_corner_r = tmp_as.mouth_corner_r(validating, :);
annotation_struct.mouth_corner_l = tmp_as.mouth_corner_l(validating, :);
annotation_struct.nose = tmp_as.nose(validating, :);
annotation_struct.N = numel(annotation_struct.names);
% save('./MAT/validating.mat', 'image');
save('./MAT/VAL.mat', 'annotation_struct', 'indices');

% 20% (the rest) of images will form testing set
testing = indices_perm; indices = testing;
% save testing set
% image = tmp_image(testing);
annotation_struct.bbox = tmp_as.bbox(testing, :);
annotation_struct.names = tmp_as.names(testing);
annotation_struct.eye_r = tmp_as.eye_r(testing, :);
annotation_struct.eye_l = tmp_as.eye_l(testing, :);
annotation_struct.canthus_rr = tmp_as.canthus_rr(testing, :);
annotation_struct.canthus_rl = tmp_as.canthus_rl(testing, :);
annotation_struct.canthus_lr = tmp_as.canthus_lr(testing, :);
annotation_struct.canthus_ll = tmp_as.canthus_ll(testing, :);
annotation_struct.mouth = tmp_as.mouth(testing, :);
annotation_struct.mouth_corner_r = tmp_as.mouth_corner_r(testing, :);
annotation_struct.mouth_corner_l = tmp_as.mouth_corner_l(testing, :);
annotation_struct.nose = tmp_as.nose(testing, :);
annotation_struct.N = numel(annotation_struct.names);
% save('./MAT/testing.mat', 'image');
save('./MAT/TST.mat', 'annotation_struct', 'indices');

%% Timestamp

fprintf(1,'Finished on %s\n\n', datestr(now));

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
%            3
%1'----------------------2'
%  |                    |
%  |                    |
% 1|                    |2
%  |                    |
%  |                    |
%4'----------------------3'
%            4

%% it is better to define by mapping cuboid to rectangle
GeoRule(1).viewID = 1;
GeoRule(1).surfaceID = 2;
GeoRule(1).pointID = [6 8 4 2];
GeoRule(2).viewID = 2;
GeoRule(2).surfaceID = 1;
GeoRule(2).pointID = [7 5 1 3];
GeoRule(3).viewID = 3;
GeoRule(3).surfaceID = 4;
GeoRule(3).pointID = [8 7 3 4];
GeoRule(4).viewID = 4;
GeoRule(4).surfaceID = 3;
GeoRule(4).pointID = [5 6 2 1];
GeoRule(5).viewID = 5;
GeoRule(5).surfaceID = 6;
GeoRule(5).pointID = [7 8 6 5];
GeoRule(6).viewID = 6;
GeoRule(6).surfaceID = 5;
GeoRule(6).pointID = [1 2 4 3];

CubRule(1).suf = [1 3];
CubRule(1).cont = [4 8];
CubRule(1).check = 5;

CubRule(2).suf = [4 1];
CubRule(2).cont = [2 6];
CubRule(2).check = 5;

CubRule(3).suf = [2 3];
CubRule(3).cont = [3 7];
CubRule(3).check = 5;

CubRule(4).suf = [2 4];
CubRule(4).cont = [1 5];
CubRule(4).check = 5;

CubRule(5).suf = [5 1];
CubRule(5).cont = [6 8];
CubRule(5).check = [3 4];

CubRule(6).suf = [5 2];
CubRule(6).cont = [5 7];
CubRule(6).check = [3 4];

CubRule(7).suf = [3 5];
CubRule(7).cont = [7 8];
CubRule(7).check = [1 2];

CubRule(8).suf = [4 5];
CubRule(8).cont = [5 6];
CubRule(8).check = [1 2];

%% cuboid rule version 1
CuboidRule(1).lhs.vid = 1;
CuboidRule(1).lhs.lid = 2;
CuboidRule(1).lhs.pid = [3 2];
CuboidRule(1).rhs.vid = 3;
CuboidRule(1).rhs.lid = 1;
CuboidRule(1).lhs.pid = [4 1];
CuboidRule(1).chk.condition = [2 3 0];
CuboidRule(1).chk.vid = 5;
CuboidRule(1).chk.lhslid = [3 2];
CuboidRule(1).chk.rhslid = [3 3];

CuboidRule(2).lhs.vid = 3;
CuboidRule(2).lhs.lid = 2;
CuboidRule(2).rhs.vid = 2;
CuboidRule(2).rhs.lid = 1;
CuboidRule(2).chk.condition = [3 0];
CuboidRule(2).chk.vid = 5;
CuboidRule(2).chk.lhslid = [3 3];
CuboidRule(2).chk.rhslid = [3 1];

CuboidRule(3).lhs.vid = 4;
CuboidRule(3).lhs.lid = 2;
CuboidRule(3).rhs.vid = 1;
CuboidRule(3).rhs.lid = 1;
CuboidRule(3).chk.condition = [3 0];
CuboidRule(3).chk.vid = 5;
CuboidRule(3).chk.lhslid = [3 4];
CuboidRule(3).chk.rhslid = [3 2];

CuboidRule(4).lhs.vid = 2;
CuboidRule(4).lhs.lid = 2;
CuboidRule(4).rhs.vid = 4;
CuboidRule(4).rhs.lid = 1;
CuboidRule(4).chk.condition = [3 0];
CuboidRule(4).chk.vid = 5;
CuboidRule(4).chk.lhslid = [3 1];
CuboidRule(4).chk.rhslid = [3 4];

%-------------------------------------
CuboidRule(5).lhs.vid = 1;
CuboidRule(5).lhs.lid = 3;
CuboidRule(5).rhs.vid = 5;
CuboidRule(5).rhs.lid = 2;
CuboidRule(5).chk(1).condition = [2 1];
CuboidRule(5).chk(1).vid = 4;
CuboidRule(5).chk(1).lhslid = [1 2];
CuboidRule(5).chk(1).rhslid = [4 3];
CuboidRule(5).chk(2).condition = [2 0];
CuboidRule(5).chk(2).vid = 3;
CuboidRule(5).chk(2).lhslid = [2 1];
CuboidRule(5).chk(2).rhslid = [3 3];

CuboidRule(6).lhs.vid = 2;
CuboidRule(6).lhs.lid = 3;
CuboidRule(6).rhs.vid = 5;
CuboidRule(6).rhs.lid = 1;
CuboidRule(6).chk(1).condition = [2 1];
CuboidRule(6).chk(1).vid = 4;
CuboidRule(6).chk(1).lhslid = [2 1];
CuboidRule(6).chk(1).rhslid = [3 1];
CuboidRule(6).chk(2).condition = [2 0];
CuboidRule(6).chk(2).vid = 3;
CuboidRule(6).chk(2).lhslid = [1 2];
CuboidRule(6).chk(2).rhslid = [3 3];

CuboidRule(7).lhs.vid = 3;
CuboidRule(7).lhs.lid = 3;
CuboidRule(7).rhs.vid = 5;
CuboidRule(7).rhs.lid = 3;
CuboidRule(7).chk(1).condition = [1 0];
CuboidRule(7).chk(1).vid = 1;
CuboidRule(7).chk(1).lhslid = [1 2];
CuboidRule(7).chk(1).rhslid = [2 3];
CuboidRule(7).chk(2).condition = [1 1];
CuboidRule(7).chk(2).vid = 2;
CuboidRule(7).chk(2).lhslid = [2 1];
CuboidRule(7).chk(2).rhslid = [1 3];

CuboidRule(8).lhs.vid = 4;
CuboidRule(8).lhs.lid = 3;
CuboidRule(8).rhs.vid = 5;
CuboidRule(8).rhs.lid = 4;
CuboidRule(8).chk(1).condition = [1 0];
CuboidRule(8).chk(1).vid = 1;
CuboidRule(8).chk(1).lhslid = [2 1];
CuboidRule(8).chk(1).rhslid = [2 3];
CuboidRule(8).chk(2).condition = [1 1];
CuboidRule(8).chk(2).vid = 2;
CuboidRule(8).chk(2).lhslid = [1 2];
CuboidRule(8).chk(2).rhslid = [1 3];

author: zouchuhang

repo: LayoutNet


text:
function [D1]=triSurfSetDist(F1,V1,F2,V2,distMetric)

% function [D1]=triSurfSetDist(F1,V1,F2,V2,distMetric)
% ------------------------------------------------------------------------
%
% Change log:
% 2021/07/22 KMM Removed waitbar for ray-tracing method
% ------------------------------------------------------------------------

%%

switch distMetric
    case 'dist'
        D1=minDist(V1,V2);
    case 'ray'
        [~,~,N1]=patchNormal(F1,V1);

        %Ray trace 1 onto 2
        optionStructRayTrace.tolEps = 1e-6;
        optionStructRayTrace.triSide = 0;
        optionStructRayTrace.rayType = 'ray'; % or line
        optionStructRayTrace.exclusionType = 'normal'; % or exclusive
        optionStructRayTrace.paired=0; % 1 for paired, 0 for non-paired

        [~,indIntersect,d1_all]=triSurfRayTrace(V1,N1,F2,V2,optionStructRayTrace);

        numIntersect=size(indIntersect,1);

        S=sparse(indIntersect(:,1),(1:numIntersect)',1,size(V1,1),numIntersect,numIntersect);
        d1_sparse=sparse(indIntersect(:,1),(1:numIntersect)',abs(d1_all),size(V1,1),numIntersect,numIntersect);

        D1=full(spmin(d1_sparse,[],2,'includenan',S~=0,1));

    case 'dist-ray'
        [D1d]=triSurfSetDist(F1,V1,F2,V2,'dist');
        [D1r]=triSurfSetDist(F1,V1,F2,V2,'ray');
        D1=gnanmin([D1d D1r],[],2);
    case 'near-norm'
        [~,indMin]=minDist(V1,V2);
        [~,~,Nv]=patchNormal(F2,V2);
        N=Nv(indMin,:);
        W=V1-V2(indMin,:);
        D1=abs(dot(N,W,2));
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
classdef chebguiExporterIVP < chebguiExporter
%CHEBGUIEXPORTERIVP   Export an IVP from CHEBGUI.
%   This is a an concrete implementation of the class CHEBGUIEXPORTER, which
%   exports IVPs from CHEBGUI to .m-files, to the workspace, or to a .chebgui
%   file. It is not intended to be called directly by the end user.
%
% See also CHEBGUI, CHEBGUIEXPORTER.

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

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS PROPERTIES:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties

        % The default file name when exporting to an .m-file:
        defaultFileName = 'chebivp.m';

        % Description for printing to .m files:
        description = 'an initial-value problem';

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function A = chebguiExporterIVP(varargin)
            % Do nothing!
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% STATIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = true )

        % Extract information from the CHEBGUI object to a struct
        expInfo = exportInfo(guifile)

        % Print problem description:
        printDescription(fid, expInfo)

        % Print options for solving the problem:
        printOptions(fid, expInfo)

        % Print lines for setting up the problem:
        printSetup(fid, expInfo, guifile)

        % Print the actual lines for calling the solver method:
        printSolver(fid, expInfo)

        % Print steps taken after the solver finishes:
        printPostSolver(fid, expInfo)

        function toWorkspaceSolutionOnly(handles)
        %TOWORKSPACESOLUTIONONLY   Export the solution to the workspace
        %
        % Note: this method exports fewer objects to the workspace than the
        % toWorkspace() method.

            % Obtain the variable names:
            varnames = handles.varnames;
            nv = numel(varnames);

            % Obtain the latest solution:
            sol = handles.latest.solution;

            % Export to the workspace!
            for k = 1:nv
                assignin('base', varnames{k}, sol(:,k));
                evalin('base', varnames{k});
            end
        end

        function toMat(handles)
        %TOMAT   Export from the CHEBGUI figure to a .mat-file.

            % Obtain the variable names involved in the problem.
            varnames = handles.varnames;
            for k = 1:numel(varnames);
                eval([varnames{k} ' = handles.latest.solution(:,k);']);
            end

            % Extract the norm of the updates, the CHEBOP and all options.
            N = handles.latest.chebop;  %#ok<NASGU>
            options = handles.latest.options;  %#ok<NASGU>

            % Show the user a figure for selecting where to save.
            uisave([varnames', 'N', 'options'], 'ivp');
        end

        function toWorkspace(handles)
        %TOWORKSPACE   Export from the CHEBGUI figure to the workspace.

            % Setup dialog for asking the user for variable names to be used.
            numlines = 1;
            options.Resize ='on';
            options.WindowStyle ='modal';

            % Need slightly different dialogs depending on whether we were
            % solving scalar problem or a coupled system.
            varnames = handles.varnames;
            nv = numel(varnames);
            if ( nv == 1 )
                prompt = {'Differential operator', 'Solution:', 'Options'};
            else
                prompt = ['Differential operator', varnames.', 'Options'];
            end

            name = 'Export to workspace';

            defaultAnswer = ['N', varnames', 'options'];

            % Show the user a dialog.
            answer = inputdlg(prompt, name, numlines, defaultAnswer, options);

            % Obtain the latest solution:
            sol = handles.latest.solution;

            % Export to the workspace!
            if ( ~isempty(answer) )
                assignin('base', answer{1}, handles.latest.chebop);
                for k = 1:nv
                    assignin('base', answer{k+1}, sol(:,k));
                end
                assignin('base', answer{nv+2}, handles.latest.options);
            end
        end

    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [y] = ml_gmm_pdf(xs,Priors,Mu,Sigma)
%ML_GMM_PDF, Gaussian Mixture Model function
%
%   G(x;\{pi,mu,sigma\}) = \sum_{i=1}^{K} pi_{k} * g(x;mu_k,sigma_k)
%
%   input ---------------------------------------------------------------
%
%       o xs: (D x N), set of N data points of dimension D
%
%       o Priors: (1 x K), sum(Priors) = 1
%
%       o Mu: (D x K), set of Sigmas
%
%       o Sigmas: (D x D x K)
%
%   output ----------------------------------------------------------------
%
%       o ys: (N x 1), values
%

[N,M] = size(xs);
K = size(Priors,2);
ys = zeros(M,K);

if size(Mu,1) == 1
    % Univariate

    for k=1:K
        ys(:,k) = normpdf(xs,Mu(k),sqrt(Sigma(k)));
    end
else
    % Multivariate
    for k=1:K
        ys(:,k) = ml_gaussPDF(xs,Mu(:,k),Sigma(:,:,k));
    end

end

y = sum(ys .* repmat(Priors,M,1),2);

end

author: epfl-lasa


text:
% EX_MAXWELL_EIG_SQUARE: solve Maxwell eigenproblem in the unit square.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_square.txt';

% Type of boundary conditions
problem_data.nmnn_sides   = [];
problem_data.drchlt_sides = [1 2 3 4];

% Physical parameters
problem_data.c_elec_perm = @(x, y) ones(size(x));
problem_data.c_magn_perm = @(x, y) ones(size(x));

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.degree     = [3 3];     % Degree of the bsplines
method_data.regularity = [2 2];     % Regularity of the splines
method_data.nsub       = [10 10];     % Number of subdivisions
method_data.nquad      = [4 4];     % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space, eigv, eigf] = ...
                             solve_maxwell_eig (problem_data, method_data);

% 4) POSTPROCESSING
[eigv, perm] = sort (eigv);
nzeros = numel (find (eigv < 1e-10));

fprintf ('Number of zero eigenvalues: %i \n', nzeros)
fprintf ('First nonzero eigenvalues: \n')
disp (eigv(nzeros+1:nzeros+5))

figure
sp_plot_solution (eigf(:,perm(nzeros+8)), space, geometry, [30 30])
title ('8^{th} eigenfunction')

%!demo
%! ex_maxwell_eig_square

%!test
%! problem_data.geo_name = 'geo_square.txt';
%! problem_data.nmnn_sides   = [];
%! problem_data.drchlt_sides = [1 2 3 4];
%! problem_data.c_elec_perm = @(x, y) ones(size(x));
%! problem_data.c_magn_perm = @(x, y) ones(size(x));
%! method_data.degree     = [3 3];     % Degree of the bsplines
%! method_data.regularity = [2 2];     % Regularity of the splines
%! method_data.nsub       = [10 10];     % Number of subdivisions
%! method_data.nquad      = [4 4];     % Points for the Gaussian quadrature rule
%! [geometry, msh, space, eigv, eigf] = solve_maxwell_eig (problem_data, method_data);
%! [eigv, perm] = sort (eigv);
%! nzeros = numel (find (eigv < 1e-10));
%! assert (msh.nel, 100)
%! assert (space.ndof, 312)
%! assert (nzeros, 121)
%! assert (eigv(nzeros+(1:5))/pi^2, [1.00000003326399; 1.00000003326400; 2.00000006652799; 4.00000968384168; 4.00000968384168], 2e-14)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [sys,values] = double(varargin)
% OBSOLETE. Use VALUE instead

% New syntax
switch nargout
    case 0
        double(value(varargin{:}))
    case 1
        sys = double(value(varargin{:}));
    case 2
        [sys,values] = value(varargin{:});
        sys = double(sys);
        values = double(values);
    otherwise
        error('Too many output arguments.');
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% AssignmentToIndex Convert assignment to index.
%
%   I = AssignmentToIndex(A, D) converts an assignment, A, over variables
%   with cardinality D to an index into the .val vector for a factor.
%   If A is a matrix then the function converts each row of A to an index.
%
%   See also IndexToAssignment.m
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function I = AssignmentToIndex(A, D)

D = D(:)'; % ensure that D is a row vector
if (any(size(A) == 1)),
    I = cumprod([1, D(1:end - 1)]) * (A(:) - 1) + 1;
else
    I = sum(repmat(cumprod([1, D(1:end - 1)]), size(A, 1), 1) .* (A - 1), 2) + 1;
end;

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function y =fbtrim(NoisyData,Param)

%   Author(s): Farhad Bayat,
%   Email: fbayat@ee.iust.ac.ir
%   Copyright 2000-2008

% Note:
%   In some applications it is necessary to differentiate a signal such as
%   position to obtain velocity signal. But always derivation makes noise
%   in output signal and using common filters yeld delay in output signal.
%   The "fbtrim.m" provides a heuristic soulotion for this problem. It
%   contains several tunning parameters you can adjust to get suitable response.
%   NoisyData: is the Data to be trim.

%   fbtrim(NoisyData,Param) or  fbtrim(NoisyData)
%   Param:  contain the trimming preferences described in the folowing:
%   Param= [NPF,CF,MSV,LCF,Nd];
%   Leave param for default values.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%   Filter Parameters:
%
% NPF=0.5;        % Noise Power Factor
% CF=0.6;         % Correction Factor
% MSV=0.01;       % minnimum data value
% LCF=5;          % Level Change factor
% DF=5 ;          % Delay Factor
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Filter parameters' default values:
NPF=0.5;        % Noise Power Factor
CF=0.6;         % Correction Factor
MSV=0.01;       % minnimum data value
LCF=5;          % Level Change factor
DF=5 ;          % Delay Factor

if nargin==2
    NPF=Param(1);
    CF=Param(2);
    MSV==Param(3);
    LCF==Param(4);
    DF==Param(5);
end

flag=0;
dv=5*MSV;

Trimed(1)=NoisyData(1);
Dyn_1=0;

for i=2:length(NoisyData)
    Dyn=(NoisyData(i)-Trimed(i-1));
    if abs(Dyn)>DF*abs(Dyn_1)
       flag=(flag+sign(Dyn))*(1+sign(Dyn*Dyn_1))/2;
    end
    if abs(Dyn)>(1+NPF)*abs(Dyn_1)
        if flag>=LCF || flag<=-LCF
            flag=0;
            Trimed(i)=NoisyData(i);
        elseif abs(Dyn_1)<MSV
            Trimed(i)=Trimed(i-1)+sign(Dyn)*dv;
        else
            Trimed(i)=Trimed(i-1)+CF*sign(Dyn)*abs(Dyn_1);
        end
    else
        Trimed(i)= NoisyData(i);
    end
    Dyn_1=(Trimed(i)-Trimed(i-1));
end

y=Trimed;

author: Sable

repo: mcbench-benchmarks


text:
function r = mtimes(a,b,varargin)
% r = a * b

% multiplication with -1 -> inversion
if isnumeric(a)

  if size(a,2) == size(b,1)

    r = b;
    r.a = a * b.a;
    r.b = a * b.b;
    r.c = a * b.c;
    r.d = a * b.d;
    r.i = a * b.i;

    r = r.normalize;
    r.i = r.i>0.5;

  else

    assert(all(abs(a(:))==1),'Rotations can be multiplied only by 1 or -1');
    tmp = rotation.id(size(a));
    tmp.i = (1-a)./2;
    a = tmp;
  end

elseif isa(b,'vector3d')

  % apply rotation
  r = rotate_outer(b,a);

elseif isa(b,'quaternion')

  r = mtimes@quaternion(a,b,varargin{:});

else

  r = rotate_outer(b,a);

end

end


text:
function spline_test23 ( )

%*****************************************************************************80
%
%% TEST23 tests SPLINE_OVERHAUSER_VAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author
%
%    John Burkardt
%
  ndata = 4;
  ndim = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST23\n' );
  fprintf ( 1, '  SPLINE_OVERHAUSER_VAL evaluates the\n' );
  fprintf ( 1, '    Overhauser spline.\n' );
%
%  Set the data.
%
  tdata(1) = 1.0;
  ydata(1,1) =   0.0;
  ydata(2,1) =   0.0;

  tdata(2) = 2.0;
  ydata(1,2) =   1.0;
  ydata(2,2) =   1.0;

  tdata(3) = 3.0;
  ydata(1,3) =   2.0;
  ydata(2,3) = - 1.0;

  tdata(4) = 4.0;
  ydata(1,4) =   3.0;
  ydata(2,4) =   0.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data to be interpolated:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of data values = %d\n', ndata );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       T             Y\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ndata
    fprintf ( 1, '%14f  ', tdata(i) );
    for j = 1 : ndim
      fprintf ( 1, '%14f  ',  ydata(j,i) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Now evaluate the spline all over the place.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  T, Spline value\n' );
  fprintf ( 1, '\n' );

  for i = 0 : 6 * ndata + 3

    tval = ( i ) / 6.0;
    yval = spline_overhauser_val ( ndim, ndata, tdata, ydata, tval );

    fprintf ( 1, '%14f  ', tval );
    for j = 1 : ndim
      fprintf ( 1, '%14f  ',  yval(j) );
    end
    fprintf ( 1, '\n' );

  end

  return
end


text:
num = size(AllFeature1,2);
% F1 = AllFeature1(1:512,:)';
F1 = max(AllFeature1(1:512,:)', AllFeature1(513:1024,:)') / 20;%(1:512,:)
% F1 = bsxfun(@minus, F1, min(F1,[],2));
% F1 = bsxfun(@rdivide, F1, max(F1,[],2));
% F1 = bsxfun(@minus, F1, mean(F1,1));
% F1 = bsxfun(@rdivide, F1, sqrt(sum(F1.^2,2)));
% F1 = bsxfun(@minus,F1,PCAmap.mean);
% F1 = F1 * PCAmap.M;

% F2 = AllFeature2(1:512,:)';
F2 = max(AllFeature2(1:512,:)', AllFeature2(513:1024,:)') / 20;%(1:512,:)
% F2 = bsxfun(@minus, F2, min(F2,[],2));
% F2 = bsxfun(@rdivide, F2, max(F2,[],2));
% F2 = bsxfun(@minus, F2, mean(F2,1));
% F2 = bsxfun(@rdivide, F2, sqrt(sum(F2.^2,2)));
% F2 = bsxfun(@minus,F2,PCAmap.mean);
% F2 = F2 * PCAmap.M;

% F1 = AllFeature1';
% F2 = AllFeature2';
thresh2 = zeros(num,1);
for i = 1:num
%     thresh2(i) = F1(i,:) * mapping.A * F1(i,:)' + F2(i,:) * mapping.A * F2(i,:)' - 2 * F1(i,:) * mapping.G * F2(i,:)';
    thresh2(i) = pdist2(F1(i,:),F2(i,:));
%     thresh2(i) = F1(i,:) * F2(i,:)';
end;
figure;
hist(thresh2(1:3000),500);
figure;
hist(thresh2(3001:end),500);

accuracies = zeros(10,1);
for i=1:10
    test_idx = [(i-1) * 300 + 1 : i*300, (i-1) * 300 + 3001 : i*300 + 3000];
    train_idx = 1:6000;
    train_idx(test_idx) = [];
    bestc=256;
    same_label = ones(6000,1);
    same_label(3001:6000) = 0;
    % predicted_label = predict(double(lfw_label),sparse(thresh2),model);
    cmd = [' -t 0 -h 0'];
    model = svmtrain(same_label(train_idx),thresh2(train_idx),cmd);
    % model = svmtrain(double(sim_label),thresh,cmd);
    [class, accuracy, deci] = svmpredict(same_label(test_idx),thresh2(test_idx),model);
    accuracies(i) = accuracy(1);
end;
mean(accuracies)
% cmd = [' -t 0 -h 0'];
% model = svmtrain(same_label,thresh2,cmd);
% [class, accuracy, deci] = svmpredict(same_label,thresh2,model);

author: happynear

repo: FaceVerification


text:
function [y,fs,wmode,fidx]=readaif(filename,mode,nmax,nskip)
%READAIF  Read a .AIF format sound file [Y,FS,WMODE,FIDX]=(FILENAME,MODE,NMAX,NSKIP)
%
% Input Parameters:
%
%	FILENAME gives the name of the file (with optional .AIF extension) or alternatively
%                 can be the FIDX output from a previous call to READAIF
%	MODE		specifies the following (*=default):
%
%    Scaling: 's'    Auto scale to make data peak = +-1
%             'r'    Raw unscaled data (integer values)
%             'q'    Scaled to make 0dBm0 be unity mean square
%             'p' *	Scaled to make +-1 equal full scale
%             'o'    Scale to bin centre rather than bin edge (e.g. 127 rather than 127.5 for 8 bit values)
%                     (can be combined with n+p,r,s modes)
%             'n'    Scale to negative peak rather than positive peak (e.g. 128.5 rather than 127.5 for 8 bit values)
%                     (can be combined with o+p,r,s modes)
%   File I/O: 'f'    Do not close file on exit
%             'd'    Look in data directory: voicebox('dir_data')
%
%	NMAX     maximum number of samples to read (or -1 for unlimited [default])
%	NSKIP    number of samples to skip from start of file
%               (or -1 to continue from previous read when FIDX is given instead of FILENAME [default])
%
% Output Parameters:
%
%	Y        data matrix of dimension (samples,channels)
%	FS       sample frequency in Hz
%	WMODE    mode string needed to recreate the data file
%	FIDX     Information row vector containing the element listed below.
%
%           (1)  file id
%				(2)  current position in file
%           (3)  dataoff	byte offset in file to start of data
%           (4)  nsamp	number of samples
%           (5)  nchan	number of channels
%           (6)  nbyte	bytes per data value
%           (7)  bits	number of bits of precision
%           (8)  code	Data format: 1=PCM, 2=ADPCM, 6=A-law, 7=Mu-law (always 1 for AIFF)
%           (9)  fs	sample frequency
%           (10) offset
%           (11) block size
%
%   If no output parameters are specified, header information will be printed.
%
%   For stereo data, y(:,1) is the left channel and y(:,2) the right

%  future enhancements: (1) handle LIST and CAT files as well
%                       (2) deal with offset and blocksize properly
%                       (3) handle other chunk types sensibly

%	   Copyright (C) Mike Brookes 1998
%      Version: $Id: readaif.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<1 error('Usage: [y,fs,wmode,fidx]=READAIFF(filename,mode,nmax,nskip)'); end

if nargin<2 mode='p';
else mode = [mode(:).' 'p'];
end
k=find((mode>='p') & (mode<='s'));
mno=all(mode~='o');                      % scale to input limits not output limits
sc=mode(k(1));

info=zeros(1,11);
if ischar(filename)
    if any(mode=='d')
        filename=fullfile(voicebox('dir_data'),filename);
    end
    fid=fopen(filename,'rb','b');
    if fid == -1
        fn=[filename,'.aif'];
        fid=fopen(fn,'rb','b');
        if fid ~= -1 filename=fn; end
    end
    if fid == -1 error(sprintf('Can''t open %s for input',filename)); end
    info(1)=fid;
else
    info=filename;
    fid=info(1);
end

if ~info(3)
    fseek(fid,0,-1);						% read file header
    header=fread(fid,4,'uchar');
    if header' ~= 'FORM' fclose(fid); error(sprintf('File %s does not begin with a FORM group ID',filename)); end
    filen=fread(fid,1,'ulong')-4;
    header=fread(fid,4,'uchar');
    if header' ~= 'AIFF' fclose(fid); error(sprintf('File %s does not begin with a AIFF type ID',filename)); end

    fmt=0;
    data=0;
    while filen>=8						% loop to read all chunks
        header=fread(fid,4,'char');
        len=fread(fid,1,'ulong');
        lenx=len+rem(len,2);    % round up to an even number
        filen=filen-lenx-8;
        %fprintf(1,'%s chunk, %d bytes, %d bytes remaining\n',char(header'),len,filen);
        if header' == 'COMM'					% ******* found COMM chunk *********
            fmt=1;
            info(5)=fread(fid,1,'ushort');			%number of channels
            info(4)=fread(fid,1,'ulong');			%number of samples
            info(7)=fread(fid,1,'ushort');				% bits per sample
            info(6)=ceil(info(7)/8);            % bytes per sample
            info(8)=1;
            ieeex=fread(fid,5,'ushort');               % read a 10-byte extended ieee format number
            fsign=ieeex(1)>32767;
            ieeex(1)=ieeex(1)-32768*fsign;
            if ~ieeex info(9)=0;
            elseif ieeex(1)==32767 info(9)=nan;
            else
                info(9)=(1-2*fsign)*sum(pow2(ieeex(2:5)',ieeex(1)-(16398:16:16446)));
            end
            fseek(fid,lenx-18,0);				% skip to end of header
        elseif header' == 'SSND'				% ******* found DATA chunk *********
            info(10)=fread(fid,1,'ulong');			%number of channels
            info(11)=fread(fid,1,'ulong');			%number of channels
            info(3)=ftell(fid);
            data=1;
            fseek(fid,lenx-8,0);				% skip to end of chunk
        else							% ******* found unwanted chunk *********
            fseek(fid,lenx,0);
        end
    end
    if filen~=0
        fprintf(2,'READAIF warning: Inconsistent file length - %d extra bytes\n',filen);
    end
end
fs=info(9);
if ~fmt | ~data fclose(fid); error(sprintf('File %s does not contain COMM and SSND chunks',filename)); end

if nargin<4 nskip=info(2);
elseif nskip<0 nskip=info(2);
end

ksamples=info(4)-nskip;
if nargin>2
    if nmax>=0
        ksamples=min(nmax,ksamples);
    end
elseif ~nargout
    ksamples=min(5,ksamples);
end
if ksamples>0
    info(2)=nskip+ksamples;
    pk=pow2(0.5,8*info(6))*(1+(mno/2-all(mode~='n'))/pow2(0.5,info(7)));  % use modes o and n to determine effective peak
    fseek(fid,info(3)+info(6)*info(5)*nskip,-1);
    nsamples=info(5)*ksamples;
    if info(6)<3
        if info(6)<2
            y=fread(fid,nsamples,'schar');
        else
            y=fread(fid,nsamples,'short');
        end
    else
        if info(6)<4
            y=fread(fid,3*nsamples,'uchar');
            y=reshape(y,3,nsamples);
            y=[1 256 65536]*y-pow2(fix(pow2(y(3,:),-7)),24);
        else
            y=fread(fid,nsamples,'long');
        end
    end
    if sc ~= 'r'
        if sc=='s' sf=1/max(max(abs(y)),1);
        elseif sc=='p' sf=1/pk;
        else   sf=2.03033976/pk;
        end
        y=sf*y;
    else
        y=y*pow2(1,info(7)-8*info(6)); % shift to get the data into the LSB end
    end

    if info(5)>1 y = reshape(y,info(5),ksamples).'; end
else
    y=[];
end

if all(mode~='f') fclose(fid); end

if nargout>2
    if ~mno wmode=[wmode 'o']; end
    if any(mode=='n') wmode=[wmode 'n']; end
    wmode=[sc num2str(info(7))];
    fidx=info;
elseif ~nargout
    fprintf(1,'\n%d Hz sample rate\n%d channel x %d samples x %d bit = %.3g seconds\n',info([9 5 4]), info(7),info(4)/info(9));
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [P,Gvec] = autoGen_getPoints(q1,q2,q3,q4,q5,l1,l2,l3,l4,l5,c1,c2,c3,c4,c5)
%AUTOGEN_GETPOINTS
%    [P,GVEC] = AUTOGEN_GETPOINTS(Q1,Q2,Q3,Q4,Q5,L1,L2,L3,L4,L5,C1,C2,C3,C4,C5)

%    This function was generated by the Symbolic Math Toolbox version 6.2.
%    22-Oct-2015 19:14:40

t2 = sin(q1);
t3 = cos(q1);
t4 = l1.*t3;
t5 = sin(q2);
t6 = cos(q2);
t7 = l2.*t6;
t8 = sin(q4);
t9 = l4.*t8;
t10 = cos(q4);
t11 = sin(q3);
t12 = cos(q3);
t13 = l3.*t12;
t14 = sin(q5);
t15 = cos(q5);
P = [-l1.*t2;t4;-l1.*t2-l2.*t5;t4+t7;-l1.*t2-l2.*t5-l3.*t11;t4+t7+t13;t9-l1.*t2-l2.*t5;t4+t7-l4.*t10;t9-l1.*t2-l2.*t5+l5.*t14;t4+t7-l4.*t10-l5.*t15];
if nargout > 1
    Gvec = [c1.*t2-l1.*t2;t4-c1.*t3;c2.*t5-l1.*t2-l2.*t5;t4+t7-c2.*t6;c3.*t11-l1.*t2-l2.*t5-l3.*t11;t4+t7+t13-c3.*t12;c4.*t8-l1.*t2-l2.*t5;t4+t7-c4.*t10;t9+c5.*t14-l1.*t2-l2.*t5;t4+t7-c5.*t15-l4.*t10];
end

author: MatthewPeterKelly

repo: OptimTraj


text:
% STD_ERPIMAGE - Compute ERP images and save them on disk.
%
% Usage:
%   >> std_erpimage( EEG, 'key', 'val', ...);
%
% Inputs:
%   EEG          - a loaded epoched EEG dataset structure. May be an array
%                  of such structure containing several datasets.
%
% Optional inputs:
%   'components' - [numeric vector] components of the EEG structure for which
%                  the measure will be computed {default|[] -> all}
%   'channels'   - [cell array] channels of the EEG structure for which
%                  activation ERPs will be computed {default|[] -> none}
%   'trialindices' - [cell array] indices of trials for each dataset.
%                  Default is all trials.
%   'recompute'  - ['on'|'off'] force recomputing data file even if it is
%                  already on disk.
%   'rmcomps'    - [integer array] remove artifactual components (this entry
%                  is ignored when plotting components). This entry contains
%                  the indices of the components to be removed. Default is none.
%   'interp'     - [struct] channel location structure containing electrode
%                  to interpolate ((this entry is ignored when plotting
%                  components). Default is no interpolation.
%   'fileout'    - [string] name of the file to save on disk. The default
%                  is the same name (with a different extension) as the
%                  dataset given as input.
%
% ERPimage options:
%   'concatenate' - ['on'|'off'] concatenate single trial of different
%                  subjects for plotting ERPimages ('on'). The default
%                  ('off') computes an ERPimage for each subject and then
%                  averages these ERPimages. This allows to perform
%                  statistics (the 'on' options does not allow statistics).
%   'smoothing'  - Smootmeter (number of trials). {Default: 10}
%                  ERPIMAGE equivalent: 'avewidth'
%   'nlines'     - Number of lines for ERPimage. ERPIMAGE equivalent is
%                  'decimate'. Note that this parameter must be larger than
%                  the minimum number of trials in each design cell
%                  {Default: 10}
%   'sorttype'   - Sorting event type(s) ([int vector]; []=all). See Notes below.
%                  Either a string or an integer.
%   'sortwin'    - Sorting event window [start, end] in milliseconds ([]=whole epoch)
%   'sortfield'  - Sorting field name. {default: latency}.
%   'erpimageopt'  - ERPIMAGE options, separated by commas (Ex: 'erp', 'cbar').
%                  {Default: none}. For further details see >> erpimage help
% Outputs:
%   erpimagestruct - structure containing ERPimage information that is
%                    been saved on disk.
%
%   Files are saved on disk.
%    [dataset_file].icaerpim     % component ERPimage file
% OR
%    [dataset_file].daterpim     % channel ERPimage file
%
% Author: Arnaud Delorme, SCCN & CERCO, CNRS, 2011-

% Copyright (C) 2011 Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. 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 HOLDER 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.

function allerpimage = std_erpimage( EEG, varargin)

if nargin < 1
    help std_erpimage;
    return;
end

allerpimage = [];
[opt, moreopts] = finputcheck( varargin, { ...
    'components'    'integer'     []                    [];
    'channels'      { 'cell','integer' }  { [] [] }     {};
    'trialindices' { 'integer','cell' }   []            [];
    'recompute'     'string'      { 'on','off' }        'off';
    'savefile'      'string'      { 'on','off' }        'on';
    'fileout'       'string'      []                    '';
    'rmcomps'       'cell'        []                    cell(1,length(EEG));
    'interp'        'struct'      { }                   struct([]);
    'nlines'        ''            []                    10;
    'smoothing'     ''            []                    10;
    'sorttype'       ''           {}                    '';
    'sortwin'        ''           {}                    [];
    'sortfield'      ''           {}                    'latency';
    'concatenate'   'string'      { 'on'  }             'on';
    'trialinfo'     'struct'      []                    struct([]);
    'savetrials'    'string'      { 'on','off' }        'on'; % obsolete (never used)
    'erpimageopt'   'cell'        {}                    {}}, ...
    'std_erpimage', 'ignore');
if ischar(opt), error(opt); end
if length(EEG) == 1 && isempty(opt.trialindices), opt.trialindices = { [1:EEG.trials] }; end
if isempty(opt.trialindices), opt.trialindices = cell(length(EEG)); end
if ~iscell(opt.trialindices), opt.trialindices = { opt.trialindices }; end
if isempty(opt.channels)
    if isfield(EEG,'icaweights')
        numc = size(EEG(1).icaweights,1);
    else
        error('EEG.icaweights not found');
    end
    if isempty(opt.components)
        opt.components = 1:numc;
    end
end

% filename
% --------
if isempty(opt.fileout), opt.fileout = fullfile(EEG(1).filepath, EEG(1).filename(1:end-4)); end
if ~isempty(opt.channels)
    filenameshort = [ opt.fileout '.daterpim'];
    prefix = 'chan';
    if iscell(opt.channels)
        if ~isempty(opt.interp)
            opt.indices = eeg_chaninds(opt.interp, opt.channels, 0);
        else
            opt.indices = eeg_chaninds(EEG(1), opt.channels, 0);
            for ind = 2:length(EEG)
                if ~isequal(eeg_chaninds(EEG(ind), opt.channels, 0), opt.indices)
                    error([ 'Channel information must be consistent when ' 10 'several datasets are merged for a specific design' ]);
                end
            end
        end
    else
        opt.indices = opt.channels;
    end
else
    opt.indices = opt.components;
    filenameshort = [ opt.fileout '.icaerpim'];
    prefix = 'comp';
end
filename = filenameshort;

% ERP information found in datasets
% ---------------------------------
if exist(filename) && strcmpi(opt.recompute, 'off')
    fprintf('Use existing file for ERSP: %s; check the ''recompute checkbox'' to force recomputing.\n', filenameshort);
    return;
end

allerpimage = [];
if strcmpi(opt.concatenate, 'off')
    % compute ERP images
    % ------------------
    if isempty(opt.channels)
         X = eeg_getdatact(EEG, 'component', opt.indices, 'trialindices', opt.trialindices );
    else X = eeg_getdatact(EEG, 'channel'  , opt.indices, 'trialindices', opt.trialindices, 'rmcomps', opt.rmcomps, 'interp', opt.interp);
    end
    if ~isempty(opt.sorttype)
         events = eeg_getepochevent(EEG, 'type', opt.sorttype, 'timewin', opt.sortwin, 'fieldname', opt.sortfield, 'trials', opt.trialindices);
    else events = [];
    end

    % reverse engeeneering the number of lines for ERPimage
    finallines = opt.nlines;
    if ~isempty(events)
         if all(isnan(events))
             error('Cannot sort trials for one of the dataset');
         end
         lastx  = sum(~isnan(events));
    else lastx  = size(X,3);
    end
    if lastx < finallines + floor((opt.smoothing-1)/2) + 3
        error('The default number of ERPimage lines is too large for one of the dataset');
    end
    firstx = 1;
    xwidth = opt.smoothing;
    %xadv   = lastx/finallines;
    nout   = finallines; %floor(((lastx-firstx+xadv+1)-xwidth)/xadv);
    nlines = (lastx-xwidth)/(nout-0.5)*i; % make it imaginary
    %nlines = ceil(lastx/((lastx-firstx+1-xwidth)/(nout-1)));

    if 0
        % testing conversion back and forth
        % ---------------------------------
        for lastx = 20:300
            for xwidth = 1:19
                for nlines = (xwidth+1):100

                    nout = floor(((lastx+nlines)-xwidth)/nlines);
                    realnlines = (lastx-xwidth)/(nout-0.5);
                    noutreal = floor(((lastx+realnlines)-xwidth)/realnlines);

                    if nout ~= noutreal
                        error('Wrong conversion 2');
                    end

                end
            end
        end
    end

    clear tmperpimage eventvals;
    for index = 1:size(X,1)
        [tmpX, tmpevents] = erpimage(squeeze(X(index,:,:)), events, EEG(1).times, '', opt.smoothing, nlines, 'noplot', 'on', opt.erpimageopt{:}, moreopts{:});
        if isempty(events), tmpevents = []; end
        eventvals{index}   = tmpevents;
        tmperpimage{index} = tmpX';
    end
    allerpimage.events = eventvals{1};
    for index = 1:size(X,1)
        allerpimage.([ prefix int2str(opt.indices(index)) ]) = tmperpimage{index};
    end
else
    % generate dynamic loading commands
    % ---------------------------------
    for dat = 1:length(EEG)
        filenames{dat} = fullfile(EEG(dat).filepath, EEG(dat).filename);
    end
    allerpimage.times = EEG(1).times;
    for index = 1:length(opt.indices)
        if ~isempty(opt.channels)
             com = sprintf('squeeze(eeg_getdatact(%s, ''interp'', chanlocsforinterp));', vararg2str( { filenames 'channel'  , opt.indices(index), 'rmcomps', opt.rmcomps, 'trialindices', opt.trialindices } ));
        else com = sprintf('squeeze(eeg_getdatact(%s));', vararg2str( { filenames 'component', opt.indices(index), 'trialindices', opt.trialindices } ));
        end
        allerpimage = setfield(allerpimage, [ prefix int2str(opt.indices(index)) ], com);
    end
    if ~isempty(opt.channels)
        com = sprintf('squeeze(eeg_getdatact(%s, ''interp'', chanlocsforinterp));', vararg2str( { filenames 'rmcomps', opt.rmcomps, 'trialindices', opt.trialindices } ));
        allerpimage = setfield(allerpimage, [ prefix 'all' ], com);
    end
    allerpimage = setfield(allerpimage, 'chanlocsforinterp', opt.interp);
    if ~isempty(opt.sorttype)
         events = eeg_getepochevent(EEG, 'type', opt.sorttype, 'timewin', opt.sortwin, 'fieldname', opt.sortfield, 'trials', opt.trialindices);
         %geteventcom = sprintf('eeg_getepochevent(%s);', vararg2str( { filenames 'type', opt.sorttype, 'timewin', opt.sortwin, 'fieldname', opt.sortfield } ));
    else events = [];
    end
    allerpimage = setfield(allerpimage, 'events', events);
end
allerpimage.times       = EEG(1).times;
allerpimage.parameters  = varargin;
allerpimage.datatype    = 'ERPIMAGE';
allerpimage.datafiles   = computeFullFileName( { EEG.filepath }, { EEG.filename });
allerpimage.datatrials  = opt.trialindices;
allerpimage.trialinfo   = opt.trialinfo;

% Save ERPimages in file (all components or channels)
% ----------------------------------------------
if strcmpi(opt.savefile, 'on')
    if strcmpi(prefix, 'comp')
        std_savedat(filename, allerpimage);
    else
        tmpchanlocs = EEG(1).chanlocs;
        allerpimage.labels = opt.channels;
        std_savedat(filename, allerpimage);
    end
end

% compute full file names
% -----------------------
function res = computeFullFileName(filePaths, fileNames)
for index = 1:length(fileNames)
    res{index} = fullfile(filePaths{index}, fileNames{index});
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function ellipse_monte_carlo_test01 ( )

%*****************************************************************************80
%
%% ELLIPSE_MONTE_CARLO_TEST01 uses ELLIPSE01_SAMPLE with an increasing number of points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2014
%
%  Author:
%
%    John Burkardt
%
  a = [ ...
    9.0, 1.0; ...
    1.0, 4.0 ]';
  e_test = [ ...
    0, 0; ...
    1, 0; ...
    0, 1; ...
    2, 0; ...
    1, 1; ...
    0, 2; ...
    3, 0 ]';
  r = 2.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ELLIPSE_MONTE_CARLO_TEST01\n' );
  fprintf ( 1, '  Use ELLIPSE01_SAMPLE to estimate integrals\n' );
  fprintf ( 1, '  in the ellipse x'' * A * x <= r^2.\n' );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '         N        1              X               Y  ' );
  fprintf ( 1, '             X^2               XY             Y^2             X^3\n' );
  fprintf ( 1, '\n' );

  n = 1;

  while ( n <= 65536 )

    [ x, seed ] = ellipse_sample ( n, a, r, seed );

    fprintf ( 1, '  %8d', n );

    for j = 1 : 7

      e(1:2) = e_test(1:2,j);

      value = monomial_value ( 2, n, e, x );

      result = ellipse_area1 ( a, r ) * sum ( value(1:n) ) / n;
      fprintf ( '  %14.6g', result );

    end

    fprintf ( 1, '\n' );

    n = 2 * n;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function decades_equal(hAxes,xLimits,yLimits)

  if (nargin < 2) || isempty(xLimits)
    xLimits = get(hAxes,'XLim');
  end
  if (nargin < 3) || isempty(yLimits)
    yLimits = get(hAxes,'YLim');
  end

  logScale = diff(yLimits)/diff(xLimits);
  powerScale = diff(log10(yLimits))/diff(log10(xLimits));

  set(hAxes,'Xlim',xLimits,...
            'YLim',yLimits,...
            'DataAspectRatio',[1 logScale/powerScale 1]);

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


function [filename] = FAIRwrite(T,m,filename)
[p,f,e] = fileparts(filename);
filename = fullfile(p,[f,'.jpg']);
fprintf('write data to %s\n',filename)
imwrite(fliplr(uint8(reshape(T,m)))',filename,'jpeg');
%eval(['!open ',filename])


text:
%  INTERNAL FUNCTION: Set initial conditions to the order of the solution of
%  the model (going from the Alphabetical order)
%
%

author: jmaih


text:
function opts = SRPrior(varargin)
%SRPRIOR Super-resolution image prior.
%   SRPRIOR returns a parameter structure to represent an image prior for
%   resolution using default settings.
%
%   The parameter structure consist of the following parameters:
%       - function:     Function handle to provide a regularization term
%                       for the prior.
%       - gradient:     Function handle to gradient of the regularizer.
%       - weight:       Regularization/prior weight.
%       - parameters:   Additional parameters passed to the regularizer and
%                       the gradient function.

opts = struct('function',   @huberPrior, ...        % Function handle to prior (default: Huber)
              'gradient',   @huberPrior_grad, ...   % Function handle to gradient
              'weight',     0, ...                  % Regularization weight
              'parameters', [], ...                 % Additional parameters passed to prior function/gradient
              'weightEvaluationFunction', []);

% Update with user-defined parameters
for k = 1:2:(nargin - 1)

    param = varargin{k};
    value = varargin{k+1};
    if isfield(opts, param)
        opts = setfield(opts, param, value);
    else
        error( sprintf('Invalid parameter %s', param) );
    end

end

author: thomas-koehler

repo: SupER


text:
function [bb] = tribal1(pp,ee)
%TRIBAL1 compute the circumballs associated with a 1-simplex
%triangulation embedded in R^2 or R^3.
%   [BB] = TRIBAL1(PP,EE) returns the circumscribing balls
%   associated with the edge segments in [PP,EE], such that
%   BB = [XC,YC,RC.^2].

%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 02/05/2018

    bb = pwrbal1(pp,zeros(size(pp,1),1),ee) ;

end

author: dengwirda