text:
function snail_trail( filename, callback_function )
%SNAIL_TRAIL Allows user to select AOIs repeatedly and remember/display them
%
% Example:
% snail_trail(filename, @disp)
%
% This example just prints bounds of selected areas, but you could do
% processing in the callback as well.
%
% TODO: Make this work for multi-image files.
%
% Author: Wade Schwartzkopf, NGA/R
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Setup layout
fig_hand = figure('MenuBar','none','Toolbar','none');
uip_hand = uipanel('Parent',fig_hand,'Position',[0 0 1 1],'BorderType','none');
mitm_hand = hg_mitm_viewer(uip_hand);
aoiButton = javaObjectEDT('javax.swing.JButton','Select AOI');
set(handle(aoiButton,'callbackproperties'), 'ActionPerformedCallback', ...
    @(obj,eventdata) process_aoi());

% Open file
mitm_hand.openFile(filename, true);
% Might be cleaner if we fit figure shape to image aspect ratio first like in TASER.
mitm_hand.Zoom = 'fit';
mitm_hand.PreChangeViewFcn = @() saveShapesNativeCoords();
mitm_hand.PostChangeViewFcn = @() restoreShapesLocalCoords();

% Add buttons to toolbar
main_toolbar = mitm_hand.main_toolbar;
main_toolbar.addSeparator;
main_toolbar.add(aoiButton);
main_toolbar.repaint;
main_toolbar.revalidate;

% Setup draggable rectangle for AOI selection
datasize = double([mitm_hand.Metadata{mitm_hand.Frame}.ImageData.NumCols ...
            mitm_hand.Metadata{mitm_hand.Frame}.ImageData.NumRows]);
native_coords = [datasize(1)/4 datasize(2)/4 datasize(1)/2-1 datasize(2)/2-1];
init_coords = mitm_hand.native2axescoords(native_coords(1:2));
init_coords(3:4) = native_coords(3:4)/max(mitm_hand.Zoom,1);
aoi_imrect = imrect(mitm_hand.AxesHandle, init_coords);
init_constraint = mitm_hand.native2axescoords(datasize+1);
% Force draggable rectangle to be within image, rather than entire figure/axes
setPositionConstraintFcn(aoi_imrect, makeConstrainToRectFcn('imrect', ...
    [1 init_constraint(1)], [1 init_constraint(2)]));
processed_areas = {};

    % This function shades all processed areas as reddish.
    function data_out = shade_processed_areas(data_in)
        % Doesn't currently handle muli-frame datasets, in which the masks
        % of processed area are almost certainly different for each frame.
        data_out = repmat(data_in, [1 1 3]);
        mask = zeros(size(data_in));
        for i = 1:numel(processed_areas)
            ul = mitm_hand.native2axescoords(processed_areas{i}(1:2));
            lr = ul + processed_areas{i}(3:4)/max(mitm_hand.Zoom,1) - 1;
            mask = mask | poly2mask([ul(1), lr(1), lr(1), ul(1)], ...
                [ul(2), ul(2), lr(2), lr(2)], size(data_in,1), size(data_in,2));
        end
        % Shade processed areas red
        data_out(:,:,2) = data_out(:,:,2) .* ~mask;
        data_out(:,:,3) = data_out(:,:,3) .* ~mask;
    end

    % This function uses the draggable rectangle to start AOI processing.
    function process_aoi
        % Get AOI
        aoi_screen = aoi_imrect.getPosition();
        % Convert to image coordinate space
        ulcorner = round(mitm_hand.axescoords2native(aoi_screen(1:2)));
        lrcorner = round(mitm_hand.axescoords2native(aoi_screen(1:2) + aoi_screen(3:4) - 1));
        % Return values (UL corner and size) to calling function
        callback_args{1}(1:2) = ulcorner;
        callback_args{1}(3:4) = lrcorner-ulcorner+1;
        callback_args{2} = mitm_hand.FrameSubFileIndices{mitm_hand.Frame}; % Which frame was AOI selected in?
        try
            callback_function(callback_args{1});
            processed_areas{end+1} = callback_args{1};  % Do we need to store frame as well?
            % This function never changes (only the processed_areas
            % variable it references), but setting it here, even if to the
            % same function it was before, causes the viewer to actually
            % apply it.
            mitm_hand.DataTransformFcn = @shade_processed_areas;
        end
    end

    % Rember draggable rectangle before zoom/pan
    function saveShapesNativeCoords()
        local_coords = aoi_imrect.getPosition();
        native_coords = mitm_hand.axescoords2native(local_coords(1:2));
        native_coords(3:4) = local_coords(3:4)*max(mitm_hand.Zoom,1);
    end

    % Restore draggable rectangle before zoom/pan
    % This is messier than one would hope because of the MATLAB quirks.
    function restoreShapesLocalCoords()
        oldzoom = get(zoom(fig_hand),'Enable');
        oldpan = get(pan(fig_hand),'Enable');
        delete(aoi_imrect);
        local_coords = mitm_hand.native2axescoords(native_coords(1:2));
        local_coords(3:4) = native_coords(3:4)/max(mitm_hand.Zoom,1);
        aoi_imrect = imrect(mitm_hand.AxesHandle, local_coords);
        ds = double([mitm_hand.Metadata{mitm_hand.Frame}.ImageData.NumCols ...
            mitm_hand.Metadata{mitm_hand.Frame}.ImageData.NumRows]);
        constraint = mitm_hand.native2axescoords(ds+1);
        setPositionConstraintFcn(aoi_imrect, makeConstrainToRectFcn('imrect', ...
            [1 constraint(1)], [1 constraint(2)]));
        % Restore old state which was changed when we redrew the shapes
        try % Doesn't work when initiated from zoom callback.
            zoom(fig_hand,oldzoom);
            pan(fig_hand, oldpan);
        end
    end
end

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

author: ngageoint


text:
function [out] = smoothThreshold_storage_logistic(S,Smax,r,e)
%smoothThreshold_storage_logistic Logisitic smoother for storage threshold functions.

% Copyright (C) 2018 Wouter J.M. Knoben
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

%   Smooths the transition of threshold functions of the form:
%
%   Q = { P, if S = Smax
%       { 0, if S < Smax
%
%   By transforming the equation above to Q = f(P,S,Smax,e,r):
%   Q = P * 1/ (1+exp((S-Smax+r*e*Smax)/(r*Smax)))
%
%   Inputs:
%   S       : current storage
%   Smax    : maximum storage
%   r       : [optional] smoothing parameter rho, default = 0.01
%   e       : [optional] smoothing parameter e, default 5
%
%   NOTE: this function only outputs the multiplier. This needs to be
%   applied to the proper flux utside of this function.
%
%   NOTE: can be applied for temperature thresholds as well (i.e. snow
%   modules). This simply means that S becomes T, and Smax T0.

% Check for inputs and use defaults if not provided
% NOTE: this is not very elegant, but it is more than a factor 10 faster then:
% if ~exist('r','var'); r = 0.01; end
% if ~exist('e','var'); e = 5.00; end
if nargin == 2
    r = 0.01;
    e = 5.00;
elseif nargin == 3
    e = 5.00;
end

% Calculate multiplier
Smax = max(Smax,0);   % this avoids numerical instabilities when Smax<0
if r*Smax == 0
    out = 1 ./ (1+exp((S-Smax+r*e*Smax)/(r)));
else
    out = 1 ./ (1+exp((S-Smax+r*e*Smax)/(r*Smax)));
end

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function [C] = spm_cov2corr(C)
% returns the correlation matrix given the covariance matrix
% FORMAT [R] = spm_cov2corr(C);
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_cov2corr.m 1143 2008-02-07 19:33:33Z spm $

%--------------------------------------------------------------------------
n    = length(C);
D    = sparse(1:n,1:n,sqrt(1./(diag(C) + eps)));
C    = real(D*C);
C    = real(C*D);

author: spm


text:
function varargout = plotstyle(varargin)
% VL_PLOTSTYLES  Get a plot style
%   STYLES = VL_PLOTSTYLE(INDEX) returns a cell array of PLOT() options
%   for the color/line style of index INDEX. The available styles are
%   defined the ColorOrder and LineStyleOrder properties of the
%   current axes and cycled first by color and then by line style.
%
%   If there is no current axes can be found, then VL_PLOTSTYLE() uses
%   the DefaultAxesColorOrder and DefaultLineStyleOrder root
%   properties.
%
%   If INDEXES has more than one element, then VL_PLOTSTYLE(INDEXES)
%   returns a cell array of styles of the same size of INDEXES.
%
%   Example::
%     The code fragment
%
%        set(gca,'LineStyleOrder', {'-','-.'}, 'NextPlot', 'replacechildren') ;
%        plot([1 2], [1:20 ; 1:20]') ;
%
%     and
%
%       set(gca,'NextPlot', 'replacechildren') ;
%       set(gca,'LineStyleOrder', {'-','-.'}) ;
%       for i = 1:20
%         style = vl_plotstyle(i) ;
%         plot([1 2], [i i], style{:}) ;
%         hold on ;
%       end
%
%     produce visually similar results.
%
%   See also: VL_HELP().
[varargout{1:nargout}] = vl_plotstyle(varargin{:});

author: yihui-he

repo: panorama

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


text:
clear;
clc;
close all;

%% https://www.coursera.org/learn/spacecraft-dynamics-kinematics/lecture/UFIBv/4-1-example-of-devenports-q-method

%% setuo the true attitude states;
theta_true = deg2rad([30 20 -10]);
BNtrue = angle2dcm(theta_true(1), theta_true(2), theta_true(3));

v1N = [1 0 0];
v2N = [0 0 1];

v1B_true = [BNtrue*v1N'];
v2B_true = [BNtrue*v2N'];

%% setup the measured attitude states
v1B = [0.8190 -0.5282 0.2243];
v2B = [-0.3138 -0.1584 0.9363];
% v1B = v1B / norm(v1B);
% v2B = v2B / norm(v2B);

% v1B = [0.8273 0.5541 -0.0920];
% v2B = [-0.8285 0.5522 -0.0955];
% v1N = [-0.1517 -0.9669 0.2050];
% v2N = [-0.8383 0.4494 -0.3044];

%% Devenport_Q mehold
w1 = 1;
w2 = 1;
W =[w1 w2];

vB = [v1B; v2B];
vB = vB.*W';
vN = [v1N; v2N];

%B = w1* v1B'*v1N + w2*v2B'*v2N;
B = vB'*vN;

S = B +B';
sigma = B(1,1) + B(2,2) + B(3,3);
Z = [B(2,3)-B(3,2) B(3,1)-B(1,3) B(1,2)-B(2,1)]';

K  = [sigma Z'; Z S - sigma*eye(3) ];

[V D] = eig(K);
[val index] = max(diag(D));
beta_q = V(:,index)';

DEVENPORT_Q = quat2dcm(beta_q)

%% OLAE  optimzal linear attitude estimator
W = eye(6);

d = [v1B - v1N,  v2B - v2N];

S =[ skew_symmetric(v1B + v1N); skew_symmetric(v2B + v2N)];
qBar = ((S'*W*S)^-1)  * S'*W*d';
OLAE_result = rod2dcm(qBar');
OLAE_result

author: yandld


text:
function [train_eval,test_eval] = ml_kcv_clustering_eval(X,labels,g,train,test,k,train_eval,test_eval,varargin)
%ML_CLUSTERING_EVAL
%
%  input ------------------------------------------------------------------
%
%       o   X        : (N x D), dataset.
%
%       o   labels   : (N x 1),  class labels : if classification (discrete)
%
%       o   g        : function handle, y = g(x)
%
%       o   train    : (M x 1),  set of indicies of X to use as train data.
%
%       o   test     : (P x 1),  set of indicies of X to use for testing.
%
%       o   k        : (1 x 1),  current iteration of k-fold
%
%       o   varargin : if you want to have a different set of labels to
%                       test put it here
%
%   output ----------------------------------------------------------------
%
%       o   train_eval : struct
%
%       o   test_eval  : struct

 % Evaluate classifier on train data

    M                        = ml_confusion_matrix(X(train(:),:),labels(train(:)),g);
    [A, P, R, F, FPR, TNR]   = ml_confusion_matrix_evaluation(M);

    if isempty(k)
        k = 1;
    end
    train_eval.accuracy(k)    = A;
    %     train_eval.precision(:,k) = P;
    %     train_eval.recall(:,k)    = R;
    train_eval.fmeasure(k)    = F;
    train_eval.fpr(k)         = FPR;
    train_eval.tnr(k)         = TNR;

    if (test~=0)
        % Evaluate classifier on test data
        if(length(varargin) == 0)
            M           = ml_confusion_matrix(X(test(:),:),labels(test(:)),g);
        else
            labelsTest = varargin{1};
            M           = ml_confusion_matrix(X(test(:),:),labelsTest(test(:)),g);
        end

        [A, P, R, F, FPR, TNR]   = ml_confusion_matrix_evaluation(M);

        test_eval.accuracy(k)    = A;
        %     train_eval.precision(:,k) = P;
        %     train_eval.recall(:,k)    = R;
        test_eval.fmeasure(k)    = F;
        test_eval.fpr(k)         = FPR;
        test_eval.tnr(k)         = TNR;
    end
end

author: epfl-lasa


text:
function A=imresize3d(V,scale,tsize,ntype,npad)
% This function resizes a 3D image volume to new dimensions
% Vnew = imresize3d(V,scale,nsize,ntype,npad);
%
% inputs,
%   V: The input image volume
%   scale: scaling factor, when used set tsize to [];
%   nsize: new dimensions, when used set scale to [];
%   ntype: Type of interpolation ('nearest', 'linear', or 'cubic')
%   npad: Boundary condition ('replicate', 'symmetric', 'circular', 'fill', or 'bound')
%
% outputs,
%   Vnew: The resized image volume
%
% example,
%   load('mri','D'); D=squeeze(D);
%   Dnew = imresize3d(D,[],[80 80 40],'nearest','bound');
%
% This function is written by D.Kroon University of Twente (July 2008)

% Check the inputs
if(exist('ntype', 'var') == 0), ntype='nearest'; end
if(exist('npad', 'var') == 0), npad='bound'; end
if(exist('scale', 'var')&&~isempty(scale)), tsize=round(size(V).*scale); end
if(exist('tsize', 'var')&&~isempty(tsize)),  scale=(tsize./size(V)); end

% Make transformation structure
T = makehgtform('scale',scale);
tform = maketform('affine', T);

% Specify resampler
R = makeresampler(ntype, npad);

% Anti-aliasing
if(scale<1)
   r=ceil(2.5/scale(1)); H=sinc((-r:r)*scale(1)); H=H./sum(H);
   Hx=reshape(H,[length(H) 1 1]);
   r=ceil(2.5/scale(2)); H=sinc((-r:r)*scale(2)); H=H./sum(H);
   Hy=reshape(H,[1 length(H) 1]);
   r=ceil(2.5/scale(3)); H=sinc((-r:r)*scale(3)); H=H./sum(H);
   Hz=reshape(H,[1 1 length(H)]);
   V=imfilter(imfilter(imfilter(V,Hx, 'same' ,'replicate'),Hy, 'same' ,'replicate'),Hz, 'same' ,'replicate');
end

% Resize the image volueme
A = tformarray(V, tform, R, [1 2 3], [1 2 3], tsize, [], 0);

author: yulequan

repo: HeartSeg

save_path: github-repos/MATLAB/yulequan-HeartSeg


text:
function bern = bernstein_poly_ab ( n, a, b, x )

%*****************************************************************************80
%
%% BERNSTEIN_POLY_AB evaluates at X the Bernstein polynomials based in [A,B].
%
%  Formula:
%
%    BERN(N,I)(X) = [N!/(I!*(N-I)!)] * (B-X)**(N-I) * (X-A)**I / (B-A)**N
%
%  First values:
%
%    B(0,0)(X) =   1
%
%    B(1,0)(X) = (      B-X                ) / (B-A)
%    B(1,1)(X) = (                 X-A     ) / (B-A)
%
%    B(2,0)(X) = (     (B-X)^2             ) / (B-A)^2
%    B(2,1)(X) = ( 2 * (B-X)    * (X-A)    ) / (B-A)^2
%    B(2,2)(X) = (                (X-A)^2  ) / (B-A)^2
%
%    B(3,0)(X) = (     (B-X)^3             ) / (B-A)^3
%    B(3,1)(X) = ( 3 * (B-X)^2  * (X-A)    ) / (B-A)^3
%    B(3,2)(X) = ( 3 * (B-X)    * (X-A)^2  ) / (B-A)^3
%    B(3,3)(X) = (                (X-A)^3  ) / (B-A)^3
%
%    B(4,0)(X) = (     (B-X)^4             ) / (B-A)^4
%    B(4,1)(X) = ( 4 * (B-X)^3  * (X-A)    ) / (B-A)^4
%    B(4,2)(X) = ( 6 * (B-X)^2  * (X-A)^2  ) / (B-A)^4
%    B(4,3)(X) = ( 4 * (B-X)    * (X-A)^3  ) / (B-A)^4
%    B(4,4)(X) = (                (X-A)^4  ) / (B-A)^4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the degree of the Bernstein polynomials to be used.
%    For any N, there is a set of N+1 Bernstein polynomials, each of
%    degree N, which form a basis for polynomials on [A,B].
%
%    Input, real A, B, the endpoints of the interval on which the
%    polynomials are to be based.  A and B should not be equal.
%
%    Input, real X, the point at which the polynomials are to be evaluated.
%
%    Output, real BERN(N+1,1), the values of the N+1 Bernstein polynomials at X.
%
  if ( b == a )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BERNSTEIN_POLY_AB - Fatal error!\n' );
    fprintf ( 1, '  A = B = %f\n', a );
    error ( 'BERNSTEIN_POLY_AB - Fatal error!' );
  end

  bern = zeros ( n + 1, 1 );

  if ( n == 0 )

    bern(1) = 1.0;

  elseif ( 0 < n )

    bern(1) = ( b - x ) / ( b - a );
    bern(2) = ( x - a ) / ( b - a );

    for i = 2 : n
      bern(i+1) = ( x - a ) * bern(i) / ( b - a );
      for j = i-1 : -1 : 1
        bern(j+1) = ( ( b - x ) * bern(j+1) + ( x - a ) * bern(j) ) / ( b - a );
      end
      bern(1) = ( b - x ) * bern(1) / ( b - a );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ y, m, d, f ] = jed_to_ymdf_roman ( jed )

%*****************************************************************************80
%
%% JED_TO_YMDF_ROMAN converts a JED to a Roman YMDF date.
%
%  Discussion:
%
%    The Roman calendar used here is artificial.  It is assumed to begin
%    on the Julian calendar date 1 January 753 BC, and to be simply a
%    copy of the Julian calendar, shifted by 753 years.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, integer Y, M, D, real F,
%    the YMDF date.
%
  [ y, m, d, f ] = jed_to_ymdf_julian ( jed );

  y = y_julian_to_roman ( y );

  return
end


text:
function c = minus(a,b)
	checkinputs(a,b);
	c = int64minus(int64(a),int64(b));
end

author: Sable

repo: mcbench-benchmarks


text:
function  write(freader,data,datatype)
%
%    write(freader,data,datatype) - writes data to a binary file opened with freader('filename')
%
fwrite(freader.fd,data,datatype);

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function [ ap, prec_k, recall_k ] = evaluate( actual, prediction, cutoff )
% Codes for evaluate average precision (AP), precision@k and recall@k
% given actual ranking and predicted ranking
%
% ARGS:
% actual       : the given actual ranking
% prediction   : the given predicted ranking
% cutoff       : cutoff k for precision@k and recall@k
%
% RETURN:
% ap           : the average precision (not MAP)
% prec_k       : the precision@k
% recall_k     : the recall@k

ap = ap_k(actual, prediction);
prks = pr_k(actual, prediction, cutoff);
prec_k = prks(1,:);
recall_k = prks(2,:);

function scores = ap_k(actual, prediction, k)
if nargin<3
    k=inf;
end

if length(prediction)>k
    prediction = prediction(1:k);
end
score = 0;
num_hits = 0;
for i=1:min(length(prediction), k)
    if sum(actual==prediction(i))>0 && ...
            sum(prediction(1:i-1)==prediction(i))==0
        num_hits = num_hits + 1;
        score = score + num_hits / i;
    end
end
scores = score / min(length(actual), k);

function scores = pr_k(actual, prediction, ks)
scores = zeros(2, length(ks));
for i=1:length(ks)
    k = ks(i);
    if length(prediction) > k
        pred = prediction(1:k);
    else
        pred = prediction;
    end
    num_hit = length(intersect(actual, pred));
    scores(1, i) = num_hit / length(pred); % precision
    scores(2, i) = num_hit / length(actual); % recall
end

author: graytowne

repo: caser

save_path: github-repos/MATLAB/graytowne-caser


text:
function sF = rdivide(sF1, sF2)
%
% Syntax
%   sF = sF1/sF2
%   sF = sF1/a
%   sF = a/sF1
%
% Input
%  sF1, sF2 - @S2FunHarmonic
%  a - double
%
% Output
%  sF - @S2FunHarmonic
%

if isnumeric(sF1)
  f = @(v) sF1./sF2.eval(v);
  sF = S2FunHarmonic.quadrature(f);
elseif isnumeric(sF2)
  sF = sF1.*(1./sF2);
else
  f = @(v) sF1.eval(v)./sF2.eval(v);
  sF = S2FunHarmonic.quadrature(f);
end

end

author: mtex-toolbox

repo: mtex

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


text:
%==========================================================================
% function onReceive(hObj, event)
%  Template callback function for the OpReadCallback of the socket manager.
%==========================================================================
function onReceive(hObj, event) %#ok<INUSL,IN%#ok<MSNU> USD>
global bAutoRequest

	%=========================================
	% received data and event timestamp
	%=========================================
	channel = event.channel;	%the SocketChannel that triggered the event
	data	= event.data;		%the received data (int8 array)
	time	= event.time;		%the time stamp

	numBytes	= numel(data);

	str = sprintf('read %.0f bytes on %s @%.3f(sec)', numBytes, char(channel.toString()), time);
	threadMessage('onReceive', str);

	%=========================================
	% TODO: do all of your processing here.
	%=========================================

	if ~isempty( findstr('10059',char(channel.toString())) ) %#ok<FSTR>
		disp('If you put a break point here you will see that the other channels');
		disp(' are still being processed while you wait here.');
	end

	%=========================================
	% Request more data.
	% NOTE: the socket manager runs on its own
	% thread and so you can set bAutoRequest
	% to false at anytime you like on the Matlab
	% command prompt to stop the
	% never-ending request-process loop once
	% it is enabled.
	%=========================================

	if( bAutoRequest )
		request		= echoServerTestBuffer();
		channelWrite(channel, request);
	end

	return;
end


text:
function varargout = process_cohere2( varargin )
% PROCESS_COHERE2: Compute the coherence between one signal in one file, and all the signals in another file.
%
% USAGE:  OutputFiles = process_cohere2('Run', sProcess, sInputA, sInputB)

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Coherence AxB [Deprecated]';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Connectivity';
    sProcess.Index       = 652;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Connectivity';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'results', 'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 2;
    sProcess.nMinFiles   = 1;
    sProcess.isPaired    = 1;

    % === CONNECT INPUT
    sProcess = process_corr2('DefineConnectOptions', sProcess);
    % === REMOVE EVOKED REPONSE
    sProcess.options.removeevoked.Comment = 'Remove evoked response from each trial';
    sProcess.options.removeevoked.Type    = 'checkbox';
    sProcess.options.removeevoked.Value   = 0;
    sProcess.options.removeevoked.Group   = 'input';
    % === COHERENCE METHOD
    sProcess.options.cohmeasure.Comment = {...
        ['Magnitude-squared Coherence
' ...
        '|C|^2 = |Gxy|^2/(Gxx*Gyy)'], ...
        ['Imaginary Coherence (2019)
' ...
        'IC    = |imag(C)|'], ...
        ['Lagged Coherence (2019)
' ...
        'LC    = |imag(C)|/sqrt(1-real(C)^2)'], ...
        [' Imaginary Coherence (before 2019)
' ...
        ' IC    = imag(C)^2 / (1-real(C)^2) ']; ...
        'mscohere', 'icohere2019','lcohere2019', 'icohere'};
    sProcess.options.cohmeasure.Type    = 'radio_label';
    sProcess.options.cohmeasure.Value   = 'mscohere';
    % === Overlap
    sProcess.options.overlap.Comment = 'Overlap for PSD estimation:' ;
    sProcess.options.overlap.Type    = 'value';
    sProcess.options.overlap.Value   = {50, '%', []};
    % === MAX FREQUENCY RESOLUTION
    sProcess.options.maxfreqres.Comment = 'Maximum frequency resolution:';
    sProcess.options.maxfreqres.Type    = 'value';
    sProcess.options.maxfreqres.Value   = {2,'Hz',2};
    % === HIGHEST FREQUENCY OF INTEREST
    sProcess.options.maxfreq.Comment = 'Highest frequency of interest:';
    sProcess.options.maxfreq.Type    = 'value';
    sProcess.options.maxfreq.Value   = {60,'Hz',2};
    % === OUTPUT MODE
    sProcess.options.outputmode.Comment = {'Save individual results (one file per input file)', 'Concatenate input files before processing (one file)', 'Save average connectivity matrix (one file)'};
    sProcess.options.outputmode.Type    = 'radio';
    sProcess.options.outputmode.Value   = 1;
    sProcess.options.outputmode.Group   = 'output';
end

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

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputA, sInputB) %#ok<DEFNU>
    % Input options
    OPTIONS = process_corr2('GetConnectOptions', sProcess, sInputA, sInputB);
    if isempty(OPTIONS)
        OutputFiles = {};
        return
    end

    % Metric options
    OPTIONS.Method = 'cohere';
    OPTIONS.RemoveEvoked  = sProcess.options.removeevoked.Value;
    OPTIONS.MaxFreqRes    = sProcess.options.maxfreqres.Value{1};
    OPTIONS.MaxFreq       = sProcess.options.maxfreq.Value{1};
    OPTIONS.CohOverlap    = 0.50;  % First pre-define the overlap
    OPTIONS.pThresh       = 0.05;
    OPTIONS.CohMeasure    = sProcess.options.cohmeasure.Value;

    % Change the overlap if it is specified
    if isfield(sProcess.options, 'overlap') && isfield(sProcess.options.overlap, 'Value') && ...
       iscell(sProcess.options.overlap.Value) && ~isempty(sProcess.options.overlap.Value) && ~isempty(sProcess.options.overlap.Value{1})
       OPTIONS.CohOverlap = sProcess.options.overlap.Value{1}/100 ;
    end

    % Compute metric
    OutputFiles = bst_connectivity({sInputA.FileName}, {sInputB.FileName}, OPTIONS);
end

author: brainstorm-tools


text:
function x = spgridm(levelseq)
% SPGRIDM  Compute grid points, maximum-norm-based grid
%    X = SPGRIDM(LEVELSEQ,D)  Computes the sparse grid points for
%    the given sequence of index LEVELSEQ and problem dimension D.
%    The coordinate value of dimension i is stored in column i of
%    the matrix X. One row of matrix X represents one grid point.
%    (Internal function)

% Author : Andreas Klimke, Universitaet Stuttgart
% Version: 1.2
% Date   : January 24, 2006

% Change log:
% V1.0   : Sep 24, 2003
% V1.1   : April 21, 2004
%          Altered function header to enable dimension-adaptive
%          grids. Simplified code.
% V1.2   : January 24, 2006
%          Changed data types to operate on uint arrays

% ------------------------------------------------------------
% Sparse Grid Interpolation Toolbox
% Copyright (c) 2006 W. Andreas Klimke, Universitaet Stuttgart
% Copyright (c) 2007-2008 W. A. Klimke. All Rights Reserved.
% See LICENSE.txt for license.
% email: klimkeas@ians.uni-stuttgart.de
% web  : http://www.ians.uni-stuttgart.de/spinterp
% ------------------------------------------------------------

% Get the number of levels
nlevels = uint32(size(levelseq,1));

% Get the dimension
d = uint16(size(levelseq,2));

npoints = zeros(nlevels,1,'uint32');

% Initialize sp with the total number of grid points of the
% level.
% Compute number of points
totalpoints = uint32(0);
for k = 1:nlevels;
	ntemp = uint32(1);
	for l = 1:d
		lev = levelseq(k,l);
		if lev == 0
			ntemp = ntemp * 3;
		else
			ntemp = ntemp * 2^uint32(lev);
		end
	end
	npoints(k) = ntemp;
	totalpoints = totalpoints + ntemp;
end

% index contains the index of the resulting array containing all
% subdomains of the level.
index = uint32(1);

x = zeros(totalpoints,d);

for kl = 1:nlevels
	level = double(levelseq(kl,:));
	for i = 1:d
		% compute the points, scaled to [0,1]
		if level(i) == 0
			c = [0; 0.5; 1];
		else
			c = (((1:2:(2^level(i).*2))).*2^(-1-level(i)))';
		end

		% Compute the number of grid points per dimension, store it
		% in repvec. The funny (level == 0) statement makes sure that
		% each level 0 dimension is counted as three.
		repvec = [2.^level+(level == 0).*2]';
		npoints = prod(repvec);
		repvec(i) = 1;
		c = repmat(shiftdim(c, 1-double(i)), repvec);
		x(index:index+npoints-1,i) = c(:);
	end
	index = index + npoints;
end


text:
function result = neg_sampling_objective(head_embedding, tail_embedding, head, tail, weights, a, b, negative_sample_rate, same_embedding)
%NEG_SAMPLING_OBJECTIVE Given a low-dimensional embedding and weights of
% 1-simplices of a high-dimensional simplicial complex, compute the
% associated negative sampling objective. This is the quantity that the
% original Python UMAP implementation is attempting to minimize by using
% the negative sampling technique during stochastic gradient descent (SGD).
% Note that this value is distinct from cross entropy as presented in the
% original UMAP paper.
% This calculation uses the modified smooth formula Phi for
% low-dimensional weight that is used in SGD.
% Because the computation time is O(n^2), we recommend only using this when
% n1*n2 <= 1e8.
%
% result = neg_sampling_objective(head_embedding, tail_embedding, head, tail, weights, a, b, negative_sample_rate, same_embedding)
%
% Parameters
% ----------
% n_neighbors: double (optional, default 15)
%     The size of local neighborhood (in terms of number of neighboring
%     sample points) used for manifold approximation. Larger values result
%     in more global views of the manifold, while smaller values result in
%     more local data being preserved. In general values should be in the
%     range 2 to 100.
%
% n_components: integer (optional, default 2)
%     The dimension of the space to embed into. This defaults to 2 to
%     provide easy visualization, but can reasonably be set to any integer
%     value in the range 2 to 100.
%
% metric: string or function (optional, default 'euclidean')
%     The metric to use to compute distances in high dimensional space. If
%     a string is passed, it must match a valid predefined metric. For now,
%     valid string metrics include:
%         * euclidean (or l2)
%         * manhattan (or l1)
%         * chebyshev (or linf)
%         * correlation
%         * cosine
%         * hamming
%         * jaccard
%         * mahalanobis
%         * minkowski
%         * seuclidean
%
% n_epochs: integer (optional)
%     The number of training epochs to be used in optimizing the low
%     dimensional embedding. Larger values result in more accurate
%     embeddings. If 0, a value will be selected based on the size of the
%     input dataset (200 for large datasets, 500 for small).
%
% learning_rate: double (optional, default 1)
%     The initial learning rate for the embedding optimization.
%
% init: string (optional, default 'spectral')
%     How to initialize the low dimensional embedding. Options are:
%         * 'spectral': use a spectral embedding of the fuzzy 1-skeleton
%         * 'random': assign initial embedding positions at random.
%         * An array of initial embedding positions.
%
% min_dist: double (optional, default 0.1)
%     The effective minimum distance between embedded points. Smaller
%     values will result in a more clustered/clumped embedding where nearby
%     points on the manifold are drawn closer together, while larger values
%     will result on a more even dispersal of points. The value should be
%     set relative to the "spread" value, which determines the scale at
%     which embedded points will be spread out.
%
% spread: double (optional, default 1)
%     The effective scale of embedded points. In combination with
%     "min_dist" this determines how clustered/clumped the embedded points
%     are.
%
% set_op_mix_ratio: double (optional, default 1)
%     Interpolate between (fuzzy) union and intersection as the set
%     operation used to combine local fuzzy simplicial sets to obtain a
%     global fuzzy simplicial sets. Both fuzzy set operations use the
%     product t-norm. The value of this parameter should be between 0 and
%     1; a value of 1 will use a pure fuzzy union, while 0 will use a pure
%     fuzzy intersection.
%
% local_connectivity: integer (optional, default 1)
%     The local connectivity required -- i.e. the number of nearest
%     neighbors that should be assumed to be connected at a local level.
%     The higher this value the more connected the manifold becomes
%     locally. In practice this should be not more than the local intrinsic
%     dimension of the manifold.
%
% repulsion_strength: double (optional, default 1)
%     Weighting applied to negative samples in low dimensional embedding
%     optimization. Values higher than one will result in greater weight
%     being given to negative samples.
%
% negative_sample_rate: integer (optional, default 5)
%     The number of negative samples to select per positive sample in the
%     optimization process. Increasing this value will result in greater
%     repulsive force being applied, greater optimization cost, but
%     slightly more accuracy.
%
% transform_queue_size: double (optional, default 4)
%     For transform operations (embedding new points using a trained model)
%     this will control how aggressively to search for nearest neighbors.
%     Larger values will result in slower performance but more accurate
%     nearest neighbor evaluation.
%
% a: double
%     Parameter of differentiable approximation of right adjoint functor.
%
% b: double
%     Parameter of differentiable approximation of right adjoint functor.ad
%
% Returns
% -------
% result: double
%     The value of the negative sampling objective.
%
% See also: CROSS_ENTROPY
%
%   AUTHORSHIP
%   Math Lead & Primary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Secondary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%
    n1 = size(head_embedding, 1);
    n2 = size(tail_embedding, 1);
    if n1*n2 > 1e8
        error('HALTED: MATLAB usually freezes for embeddings this large.');
    end

    n_other_points = n2;
    if same_embedding
        n_other_points = n_other_points - 1;
    end

    full_dists = pdist2(head_embedding, tail_embedding);

    Phi = ones(size(full_dists))./(1 + a*(full_dists.^(2*b)));
    if same_embedding
        Phi = Phi - diag(diag(Phi));
    end
    gap_part = sum(log(1 - Phi), 2);

    Phi_summands = weights.*(log(Phi(sub2ind(size(Phi), head, tail))) + negative_sample_rate/n_other_points*gap_part(head));

    result = -sum(Phi_summands);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function g = gauss(N)

% gauss(N)
%
% returns N normally distributed random numbers
% reference : Num. Recipes, Chapter 7.2 Normal Deviates

rsq  = [];

M = ceil(N/2);
while length(rsq) < M	% make shure we really have at least M values
	v = 2*rand(ceil(M*1.33),2)-1;		% produce more random numbers
	rsq = v(:,1).*v(:,1)+v(:,2).*v(:,2);
	ind = find((rsq >=1) | (rsq == 0));
	rsq(ind) = [];						% because we want to remove some
end

v(ind,:) = [];
v = v(1:M,:);
rsq = rsq(1:M);
fac = sqrt(-2 * log(rsq) ./ rsq);
g = [v(:,1) .* fac ; v(:,2) .* fac];
g = g(1:N);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [Izerosmooth]= filtering(I1)
mask= [0 0 1 0 0;
       0 1 2 1 0;
       1 2 -16 2 1;
       0 1 2 1 0;
       0 0 1 0 0];

   Izerosmooth= conv2(I1,mask);


text:
%% Theory of Misorientations
%
%%
% Misorientation describe the relative orientation of two crystal with
% respect to each other. Those crystal may be of the same phase or of
% different phases. Misorientation are used to describe

%% Grain Exchange Symmetry

%%
% Misorientation describes the relative orientation of two grains with
% respect to each other. Important concepts are twinnings and
% CSL (coincidence site lattice) misorientations. To illustrate this
% concept at a practical example let us first import some Magnesium EBSD
% data.

mtexdata twins silent

% use only proper symmetry operations
ebsd('M').CS = ebsd('M').CS.properGroup;

% compute grains
grains = calcGrains(ebsd('indexed'),'threshold',5*degree);
CS = grains.CS; % extract crystal symmetry

%%
% Next we plot the grains together with their mean orientation and
% highlight grain 74 and grain 85

plot(grains,grains.meanOrientation,'micronbar','off')

hold on
plot(grains([74,85]).boundary,'edgecolor','w','linewidth',2)
hold off

text(grains([74,85]),{'1','2'})

%%
% After extracting the mean orientation of grain 74 and 85

ori1 = grains(74).meanOrientation;
ori2 = grains(85).meanOrientation;

%%
% we may compute the misorientation angle between both orientations by

angle(ori1, ori2) ./ degree

%%
% Note that the misorientation angle is computed by default modulo crystal
% symmetry, i.e., the angle is always the smallest angles between all
% possible pairs of symmetrically equivalent orientations. In our example
% this means that symmetrisation of one orientation has no impact on the
% angle

angle(ori1, ori2.symmetrise) ./ degree

%%
% The misorientation angle neglecting crystal symmetry can be computed by

angle(ori1, ori2.symmetrise,'noSymmetry')./ degree

%%
% We see that the smallest angle indeed coincides with the angle computed
% before.

%% Misorientations
% Remember that both orientations ori1 and ori2 map crystal coordinates
% onto specimen coordinates. Hence, the product of an inverse orientation
% with another orientation transfers crystal coordinates from one crystal
% reference frame into crystal coordinates with respect to another crystal
% reference frame. This transformation is called misorientation

mori = inv(ori1) * ori2

%%
% In the present case the misorientation describes the coordinate transform
% from the reference frame of grain 85 into the reference frame of crystal
% 74. Take as an example the plane {11-20} with respect to the grain 85.
% Then the plane in grain 74 which alignes parallel to this plane can be
% computed by

round(mori * Miller(1,1,-2,0,CS))

%%
% Conversely, the inverse of mori is the coordinate transform from crystal
% 74 to grain 85.

round(inv(mori) * Miller(2,-1,-1,0,CS))

%% Coincident lattice planes
% The coincidence between major lattice planes may suggest that the
% misorientation is a twinning misorientation. Lets analyse whether there
% are some more alignments between major lattice planes.

%m = Miller({1,-1,0,0},{1,1,-2,0},{1,-1,0,1},{0,0,0,1},CS);
m = Miller({1,-1,0,0},{1,1,-2,0},{-1,0,1,1},{0,0,0,1},CS);

% cycle through all major lattice planes
close all
for im = 1:length(m)
  % plot the lattice planes of grains 85 with respect to the
  % reference frame of grain 74
  plot(mori * m(im).symmetrise,'MarkerSize',10,...
    'DisplayName',char(m(im)),'figSize','large','noLabel','upper')
  hold all
end
hold off

% mark the corresponding lattice planes in the twin
mm = round(unique(mori*m.symmetrise,'noSymmetry'),'maxHKL',6);
annotate(mm,'labeled','MarkerSize',5,'figSize','large','textAboveMarker')

% show legend
legend({},'location','SouthEast','FontSize',13);

%%
% we observe an almost perfect match for the lattice planes {11-20} to
% {-2110} and {1-101} to {-1101} and good coincidences for the lattice
% plane {1-100} to {0001} and {0001} to {0-661}. Lets compute the angles
% explicitly

angle(mori * Miller(1,1,-2,0,CS),Miller(2,-1,-1,0,CS)) / degree
angle(mori * Miller(1,0,-1,-1,CS),Miller(1,-1,0,1,CS)) / degree
angle(mori * Miller(0,0,0,1,CS) ,Miller(1,0,-1,0,CS),'noSymmetry') / degree
angle(mori * Miller(1,1,-2,2,CS),Miller(1,0,-1,0,CS)) / degree
angle(mori * Miller(1,0,-1,0,CS),Miller(1,1,-2,2,CS)) / degree

%% Twinning misorientations
% Lets define a misorientation that makes a perfect fit between the {11-20}
% lattice planes and between the {10-11} lattice planes

mori = orientation.map(Miller(1,1,-2,0,CS),Miller(2,-1,-1,0,CS),...
  Miller(-1,0,1,1,CS),Miller(-1,1,0,1,CS))

% the rotational axis
round(mori.axis)

% the rotational angle
mori.angle / degree

%%
% and plot the same figure as before with the exact twinning
% misorientation.

% cycle through all major lattice planes
close all
for im = 1:length(m)
  % plot the lattice planes of grains 85 with respect to the
  % reference frame of grain 74
  plot(mori * m(im).symmetrise,'MarkerSize',10,...
    'DisplayName',char(m(im)),'figSize','large','noLabel','upper')
  hold all
end
hold off

% mark the corresponding lattice planes in the twin
mm = round(unique(mori*m.symmetrise,'noSymmetry'),'maxHKL',6);
annotate(mm,'labeled','MarkerSize',5,'figSize','large')

% show legend
legend({},'location','NorthWest','FontSize',13);

%% Highlight twinning boundaries
% It turns out that in the previous EBSD map many grain boundaries have a
% misorientation close to the twinning misorientation we just defined. Lets
% Lets highlight those twinning boundaries

% consider only Magnesium to Magnesium grain boundaries
gB = grains.boundary('Mag','Mag');

% check for small deviation from the twinning misorientation
isTwinning = angle(gB.misorientation,mori) < 5*degree;

% plot the grains and highlight the twinning boundaries
plot(grains,grains.meanOrientation,'micronbar','off')
hold on
plot(gB(isTwinning),'edgecolor','w','linewidth',2)
hold off

%%
% From this picture we see that large fraction of grain boudaries are
% twinning boundaries. To make this observation more evident we may plot
% the boundary misorientation angle distribution function. This is simply
% the angle distribution of all boundary misorientations and can be
% displayed with

close all
plotAngleDistribution(gB.misorientation)

%%
% From this we observe that we have about 50 percent twinning boundaries.
% Analogously we may also plot the axis distribution

plotAxisDistribution(gB.misorientation,'contour')

%%
% which emphasises a strong portion of rotations about the (-12-10) axis.

%% Phase transitions
% Misorientations may not only be defined between crystal frames of the
% same phase. Lets consider the phases Magnetite and Hematite.

CS_Mag = loadCIF('Magnetite')
CS_Hem = loadCIF('Hematite')

%%
% The phase transition from Magnetite to Hematite is described in
% literature by {111}_m parallel {0001}_h and {-101}_m parallel {10-10}_h
% The corresponding misorientation is defined in MTEX by

Mag2Hem = orientation.map(...
  Miller(1,1,1,CS_Mag),Miller(0,0,0,1,CS_Hem),...
  Miller(-1,0,1,CS_Mag),Miller(1,0,-1,0,CS_Hem))

%%
% Assume a Magnetite grain with orientation

ori_Mag = orientation.byEuler(0,0,0,CS_Mag)

%%
% Then we can compute all variants of the phase transition by

symmetrise(ori_Mag) * inv(Mag2Hem)

%%
% and the corresponding pole figures by

plotPDF(symmetrise(ori_Mag) * inv(Mag2Hem),...
  Miller({1,0,-1,0},{1,1,-2,0},{0,0,0,1},CS_Hem))


text:
% AUSM (Liou-Steffen) scheme for one-dimensional Euler equations

% Copyright 2001 P. Wesseling
% This program and its subprograms may be freely used, modified and distributed
% under the GNU General Public License: http://www.gnu.org/copyleft/gpl.html

% Theory in Section 10.5 of:

% 	P. Wesseling: Principles of Computational Fluid Dynamics
% 	Springer, Heidelberg, 2000 ISBN 3-5453-0. XII, 642 pp.
% See http://ta.twi.tudelft.nl/nw/users/wesseling/cfdbook.html

% This program makes Figs. 10.17, 10.20 in the book

% Function called: f
% Programs called: problem_specification, Riemann

clear all

global  PRL  CRL MACHLEFT  gamma  pleft  pright  rholeft  rhoright  uleft...
	uright  tend  lambda		% lambda = dt/dx

		% .....................Input............................
gamma = 1.4; 	% Ratio of specific heats
J = 48;		% Number of grid cells
bouncon = 0;	% bouncon chooses outflow boundary conditions
		% = 0: Nothing happens: infinite domain
		% = 1: Solid wall at x = 1 with direct prescription of uwall
		% = 2: Solid wall at x = 1 with reflection b.c.
		% ....................End of input........................

gammab = 1/(gamma - 1); gam1 = gamma-1; gamgam = gamma^gamma;
problem_specification

h = 1/J;  				% Cell size
dt = lambda*h;				% Time step
n = floor(tend/dt);			% Number of time-steps

% 		Definition of grid numbering
%       x=0    					 x=1
% grid   |---o---|---o---|---o---  ...  --|---o---|
%            1   1   2   2   3           J-1  J

xcenter = h*[1:J] - h/2;		% Location of cell centers

press = zeros(size(xcenter));		% Preallocation of pressure,
rhoold = press; uold = press;		%       density and velocity
rhonew = press; mnew = press;		%       momentum
totenew = press; enthalpy = press;	%	total energy and enthalpy

for j = 1:length(xcenter)		% Initial conditions
  if xcenter(j) < 0.5, press(j) = pleft; rhoold(j) = rholeft;  uold(j) = uleft;
  else,  	     press(j) = pright;  rhoold(j) = rhoright; uold(j) = uright;
  end
end

	% Initialization of cell center variables
totenold = rhoold.*(0.5*uold.*uold + gammab*press./rhoold); % Total energy rho*E
totenleft = totenold(1); totenright = totenold(J);
mold = rhoold.*uold;					    % Momentum m
c = sqrt(gamma*press./rhoold);				    % Sound speed
mach = uold./c;						    % Mach number
enthalpy =  0.5*uold.*uold + gammab*c.^2;		    % Enthalpy

machplus = mach; machminus = mach;		% Preallocation of
presplus = mach; presminus = mach;		%     	split fluxes
flux1 = zeros(J-1,1); flux2 = flux1;		% 	and Liou-Steffen
flux3 = flux1; machhalf = flux1;		%	fluxes
m1 = flux1; m2 = flux1;

t = 0;
for i = 1:n,  t = t + dt;
  for j = 1:J
    if mach(j) > 1
      machplus(j) = mach(j);  machminus(j) = 0;
      presplus(j) = press(j); presminus(j) = 0;
    elseif mach(j) < -1
      machplus(j) = 0; machminus(j) = mach(j);
      presplus(j) = 0; presminus(j) = press(j);
    else
      machplus(j)  =  0.25*(mach(j) + 1)^2;
      machminus(j) = -0.25*(mach(j) - 1)^2;
      presplus(j)  =  0.5*press(j)*(1 + mach(j));
      presminus(j) =  0.5*press(j)*(1 - mach(j));
    end
  end

	% Liou-Steffen fluxes
  for j = 1:J-1,  machhalf(j) = machplus(j) + machminus(j+1); end
  m1 = 0.5*(machhalf + abs(machhalf));  m2 = 0.5*(machhalf - abs(machhalf));
  rhoc = rhoold.*c;
  for j = 1:J-1
    flux1(j) = m1(j)*rhoc(j) + m2(j)*rhoc(j+1);
    flux2(j) = m1(j)*rhoc(j)*uold(j) +  m2(j)*rhoc(j+1)*uold(j+1) +...
      presplus(j) + presminus(j+1);
    flux3(j) =  m1(j)*rhoc(j)*enthalpy(j)  + m2(j)*rhoc(j+1)*enthalpy(j+1);
  end

	% Update of state variables
  rhonew(1)  = rholeft;  	rhonew(J)  = rhoright;
  mnew(1)    = rholeft*uleft;	mnew(J)    = rhoright*uright;
  totenew(1) = totenleft; 	totenew(J) = totenright;
  for j = 2:J-1
    rhonew(j)  = rhoold(j)   - lambda*(flux1(j) - flux1(j-1));
    mnew(j)    = mold(j)     - lambda*(flux2(j) - flux2(j-1));
    totenew(j) = totenold(j) - lambda*(flux3(j) - flux3(j-1));
  end
  uold   = mnew./rhonew;   press = gam1*(totenew - 0.5*mnew.*uold);
  rhoold = rhonew; 	totenold = totenew; 	mold = mnew;
  c = sqrt(gamma*press./rhoold); 	mach = uold./c;
  enthalpy =  0.5*uold.*uold + gammab*c.^2;
end
entropy = log(press./rhoold.^gamma);

figure(1), clf
subplot(2,3,1),hold on,title('DENSITY','fontsize',14),plot(xcenter,rhonew,'o')
subplot(2,3,2),hold on,title('VELOCITY','fontsize',14),plot(xcenter,uold,'o')
subplot(2,3,3),hold on,title('PRESSURE','fontsize',14),plot(xcenter,press,'o')
subplot(2,3,4),hold on,title('MACHNUMBER','fontsize',14),plot(xcenter,mach,'o')
subplot(2,3,5),hold on,title('ENTROPY','fontsize',14),plot(xcenter,entropy,'o')
subplot(2,3,6),axis('off'),title('Liou-Steffen scheme','fontsize',14)

Riemann		% Plot exact solution

repo: Aero-matlab


text:
function[varargout]=vrep(varargin)
%VREP  Replicates an array along a specified dimension.
%
%   Y=VREP(X,N,DIM) replicates the array by N times along dimension
%   dimension DIM.  For instance:
%
%        VREP([1:4]',3,2)=[ [1:4]' [1:4]' [1:4]' ]
%
%   This is often useful in array algebra.
%
%   IF N and DIM are arrays of length M, then X is replicated along each of
%   the M different dimensions: N(1) times along dimensions DIM(1), etc.
%
%   [Y1,Y2,...,YP]=VREP(X1,X2,...,XP,N,DIM) also works.
%
%   See also VINDEX, DIM.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2001--2018 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
  vrep_test,return
end

n=varargin{end-1};
ndim=varargin{end};

for i=1:length(varargin)-2
   varargout{i}=varargin{i};
end

for i=1:length(varargin)-2
    for j=1:length(n)
        varargout{i}=vrep1(varargout{i},n(j),ndim(j));
    end
end

eval(to_overwrite(nargin-2))

%You would think Matlab would provide a simpler way to do this.
function[y]=vrep1(x,n,dim)

str='y=repmat(x,[';
ndx=ndims(x);
for i=1:max(ndx,dim)
    if i~=dim
        str=[str '1,'];
    else
	str=[str 'n,'];
    end
end
str=[str(1:end-1) ']);'];
eval(str);

function[]=vrep_test

ans1=vrep((1:4)',3,2);
ans2=[ (1:4)' (1:4)' (1:4)' ];
reporttest('VREP', aresame(ans1,ans2))

x1=(1:4)';x2=(1:4)';
vrep(x1,x2,3,2);
reporttest('VREP output redirect', aresame(x1,ans1) && aresame(x2,ans2))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [unionVolume,intersectVolume] = calcUnionIntersectionVol(structNumV,planC)
%function planC = calcUnionIntersectionVol(structNumV,planC)
%
%This function calculates volume of union and intersection from structNumV.
%
%APA, 01/25/08
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

global stateS

if ~exist('planC')
    global planC
end

indexS = planC{end};

scanNum = getStructureAssociatedScan(structNumV(1),planC);
structureC{1} = getUniformStr(structNumV(1));
%Get z-coordinates of scan
[scanXv,scanYv,scanZv] = getUniformScanXYZVals(planC{indexS.scan}(scanNum));
[numRows,numCols,numSlcs] = size(structureC{1});
zeroSlcM = zeros(numRows,numCols,'uint8');

for strInd = 2:length(structNumV)

    structNum = structNumV(strInd);

    %return if structNum is already associated to scanNum
    assocScanNum = getAssociatedScan(planC{indexS.structures}(structNum).assocScanUID,planC);

    if scanNum ~= assocScanNum
        struct3M = [];

        %Get transformation matrix
        if ~isfield(planC{indexS.scan}(scanNum),'transM') | isempty(planC{indexS.scan}(scanNum).transM)
            transMnew = eye(4);
        else
            transMnew = planC{indexS.scan}(scanNum).transM;
        end
        if ~isfield(planC{indexS.scan}(assocScanNum),'transM') | isempty(planC{indexS.scan}(assocScanNum).transM)
            transMold = eye(4);
        else
            transMold = planC{indexS.scan}(assocScanNum).transM;
        end
        transM = inv(transMnew)*transMold;

        [jnk, relStructNum] = getStructureAssociatedScan(structNum, planC);

        h = waitbar(0,'Interpolating scan...');

        for i = 1:length(scanZv)
            coord = scanZv(i);

            %structUID   = planC{indexS.structures}(structNum).strUID;
            %contourS    = calllib('libMeshContour','getContours',structUID,single(pointOnPlane),single(planeNormal),single([0 1 0]),single([1 0 0]));

            [slc, sliceXVals, sliceYVals] = getStructureSlice(assocScanNum, 3, coord, transM, planC);
            if (length(sliceXVals) ~= numCols) || (length(sliceYVals) ~= numRows)
                error('Associated scans have different size')
            end
            oneStructM  = double(bitget(slc, relStructNum));
            if ~isempty(oneStructM)
                struct3M(:,:,i) = oneStructM';
            else
                struct3M(:,:,i) = zeroSlcM;
            end
            waitbar(i/length(scanZv),h)

        end
        close(h)
        structureC{strInd} = struct3M;

    else % structure associated to same base scan

        structureC{strInd} = getUniformStr(structNum, planC);

    end

end

combineStr3M = structureC{1};
for i=2:length(structureC)
    combineStr3M = combineStr3M + structureC{i};
end
intersectStr3M = combineStr3M==length(structureC);
unionStr3M = combineStr3M>0;

dx = abs(scanXv(1)-scanXv(2));
dy = abs(scanYv(1)-scanYv(2));
dz = abs(scanZv(1)-scanZv(2));

unionVolume = length(find(unionStr3M)) * dx*dy*dz;
intersectVolume = length(find(intersectStr3M)) * dx*dy*dz;

nameStr = '';
for i=1:length(structNumV)
    nameStr = [nameStr , ', ', planC{indexS.structures}(structNumV(i)).structureName];
end
nameStr(1) = [];
nameStr = ['Structures:', nameStr];

CERRStatusString('----------------------------------------------------------','console')
CERRStatusString(nameStr,'console')
CERRStatusString(['Volume of Union = ',num2str(unionVolume),' cc'],'console')
CERRStatusString(['Volume of Intersection = ',num2str(intersectVolume), ' cc'],'console')
CERRStatusString(['Intersection/Union = ',num2str(intersectVolume/unionVolume)],'console')
CERRStatusString('----------------------------------------------------------','console')
CERRStatusString(['Intrsc = ',num2str(intersectVolume), ',  Union = ',num2str(unionVolume), ' cc'],'gui')

return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% demo operant learning with dynamic learning rate
% This demo simulates a "volatile learner" experiencing feedback to her
% chosen actions. The pseudo learning rate then changes as a function of
% (estimated) environmental volatility. We then fit the behavioural
% response with a stochastic variant of the Q-learning model, which allows
% estimating the time series of learning rates. We then perform a Volterra
% decomposition of hidden states (which include the learning rate).
% Finally, we fit an "augmented" Q-learning model, which predicts an
% transient acceleration of learning rate following changesin the winning
% action.

close all
clear variables

% simulate VB volatile learner in a 2-armed bandit task
fb.inH.u0 = repmat([ones(1,50),zeros(1,50)],1,4); % feedbacks with reversals
nt = size(fb.inH.u0,2)+1;
fb.h_fname = @h_truefalse;
fb.indy = 1;
fb.indfb = 2;
theta = [0;-2;0];
phi = [3;0]; % inverse temperature & bias
x0 = repmat([0;0;0;0;0],2,1);
inF.lev2 = 1; % 3rd level (volatility learning)
inF.kaub = 1.4;
inF.thub = 1;
inF.rf = 1;
inG.respmod = 'taylor';
options.sources.type = 1 ; % binomial observation;
options.inF = inF;
options.inG = inG;
options.skipf = zeros(1,nt);
options.skipf(1) = 1;
[y,x,x0,eta,e,u] = VBA_simulate (nt,@f_OpLearn,@g_VBvolatile0,theta,phi,zeros(2,nt),Inf,Inf,options,x0,fb);
% plot simulated behaviour
hf = figure('color',[1 1 1 ],'name','simulated choices');
ha = axes('parent',hf);
plot(ha,y-e,'r')
hold(ha,'on')
plot(ha,y,'kx')
legend(ha,{'p(y=1|theta,phi,m)','binomial data samples'})
VBA_getSubplots ();
dummy.options = options;
[ha,hf] = unwrapVBvolatileOTO(struct('muX',x,'muTheta',theta),dummy);
set(hf,'name','simulated volatile VB learner')

% dummy VB inversion (with ideal priors) of volatile learner
d00 = struct('n',2*5,'n_theta',3,'n_phi',2);
priors = [];
priors.muPhi = phi;
priors.muTheta = theta;
priors.muX0 = x0;
priors.SigmaPhi = 0*eye(d00.n_phi);
priors.SigmaTheta = 0*eye(d00.n_theta);
priors.SigmaX0 = 0*eye(d00.n);
priors.a_alpha = Inf;
priors.b_alpha = 0;
opt00 = options;
opt00.priors = priors;
[p00,o00] = VBA_NLStateSpaceModel(y,u,@f_OpLearn,@g_VBvolatile0,d00,opt00);

%----
% Now set up Q-learner model
u(fb.indfb,u(fb.indfb,:)==0)=-1; % mean-centre feedback for value learning
d0 = struct('n',3,'n_theta',0,'n_phi',1);
priors = [];
priors.a_alpha = 1;
priors.b_alpha = 1;
tmp = 1e2*eye(3);
tmp(3,3) = 1e0;
for t=1:nt
    priors.iQx{t} = tmp;
end
opt0 = [];
opt0.backwardLag = 32;
opt0.priors = priors;
opt0.sources.type = 1;
opt0.verbose = 1;
opt0.MaxIter = 3;
opt0.kernelSize = 32;
opt0.detrendU = 4;
[p0,o0] = VBA_NLStateSpaceModel(y,u,@f_Qlearn_dynLR,@g_softmax,d0,opt0);

% check relation between identified learning rate and volatility of VB-learner:
it = 1:400;
X = [VBA_vec(p0.muX(3,:)),ones(nt,1)];
Y = VBA_vec(sum(x([4,9],:),1));
[pv,stat,df,all] = GLM_contrast(X,Y,[1;0],'F',1,{'learning rate','CST'},{'volatility'});

% perform posterior Volterra decomposition
uu = zeros(3,nt);
uu(1,:) = 2*u(fb.indy,:)-1; % previous own action
for t=1:size(u,2)
    % uu(2,:) = winning action
    if u(fb.indfb,t)==1
        uu(2,t) = uu(1,t);
    else
        uu(2,t) = -uu(1,t);
    end
    % uu(3,:) = winning action stability
    try
        if ~isequal(uu(2,t),uu(2,t-1))
            uu(3,t) = 1;
        else
            uu(3,t) = 0;
        end
    end
end
o0_v = o0;
o0_v.u = uu;
[kernels] = VBA_getVolterraKernels(p0,o0_v);
hf = figure('color',[1 1 1],'name','Volterra decomposition');
mk = squeeze(kernels.x.m(3,:,:))';
vk = squeeze(kernels.x.v(3,:,:))';
ha = subplot(2,1,1,'parent',hf);
plotUncertainTimeSeries(mk,vk,[],ha)
legend(ha,{'agent''s chosen action','winning action','winning action stability'})
xlabel(ha,'lag')
ylabel(ha,'Volterra weight')
title(ha,'stochastic learning rate')

% Now fit agent's choice with augmented Q-learning model
u(3,:) = uu(3,:); % add in winning action stability for learning rate evolution
u(3,u(3,:)==0)=-1; % mean-centre for arbitrary steady-state of learning rate
d1 = struct('n',4,'n_theta',2,'n_phi',1);
priors = [];
priors.a_alpha = Inf;
priors.b_alpha = 0;
opt1.priors = priors;
opt1.sources.type = 1;
opt1.figName = 'augmented Q-learning model';
[p1,o1] = VBA_NLStateSpaceModel(y,u,@f_Qlearn_gammaLR,@g_softmax,d1,opt1);

% extract and plot Volterra kernels
o1_v = o1;
o1_v.u = uu;
[kernels] = VBA_getVolterraKernels(p1,o1_v);
mk = squeeze(kernels.x.m(3,:,:))';
vk = squeeze(kernels.x.v(3,:,:))';
ha = subplot(2,1,2,'parent',hf);
plotUncertainTimeSeries(mk,vk,[],ha)
legend(ha,{'agent''s chosen action','winning action','winning action stability'})
xlabel(ha,'lag')
ylabel(ha,'Volterra weight')
title(ha,'augmented learning rate')

author: MBB-team

repo: VBA-toolbox


text:
function L1 = calc_Laplacian_for_H1(B1)
B = B1 - 1;
L = calc_Laplacian_for_H(B);
L1 = zeros(B1,B1);
L1(2:end,2:end) = L;

function L = calc_Laplacian_for_H(B)
W = diag(ones(1,B-1),-1) + diag(ones(1,B-1),1);
D = diag(sum(W,2));
L = D - W;

author: zhouyuanzxcv

repo: Hyperspectral


text:
function [D] = spm_eeg_inv_Mesh2Voxels(varargin)
% Convert a mesh representation of M/EEG power into a smoothed image
% FORMAT [D] = spm_eeg_inv_Mesh2Voxels(D,[val])
% Input:
% D        - MEEG object or filename of M/EEG mat-file (optional)
%
%     D.inv{val}.contrast.display:   display image at the end {true, [false]}
%     D.inv{val}.contrast.space:     native [0] or MNI {1} output space
%     D.inv{val}.contrast.format:    output file format {['image'], 'mesh'}
%     D.inv{val}.contrast.smoothing: # iterations for cortical smoothing
%
% Output:
% D        - MEEG object containing the new image filenames in fields:
%
%     D.inv{val}.contrast.fname
%__________________________________________________________________________
%
% Non-linear interpolation of a Mesh contrast into MNI Voxel space
% This routine is used to produce a 3D image canonical sMRI
% space (in voxel coordinates) from a cortical mesh (3D surface).
% This yields a NIfTI image of the summary statistics of the cortical
% activity for the effect of interest. This image can then enter the
% classical SPM routines for statistical testing.
% The [non-negative] mean square contrast is smoothed both on the mesh
% (using a graph Laplacian) and then in voxel-space using a conventional
% Gaussian filter.
%__________________________________________________________________________
% Copyright (C) 2007-2017 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_eeg_inv_Mesh2Voxels.m 7094 2017-06-06 11:14:10Z guillaume $

SVNrev = '$Rev: 7094 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);

%-Parse input arguments
%--------------------------------------------------------------------------
[D,val]      = spm_eeg_inv_check(varargin{:});

%-Get options
%--------------------------------------------------------------------------
try, space   = D.inv{val}.contrast.space;     catch, space  = 1; end
try, fmt     = D.inv{val}.contrast.format;    catch, fmt    = 'image'; end
try, smooth  = D.inv{val}.contrast.smoothing; catch, smooth = 8; end
try, Disp    = D.inv{val}.contrast.display;   catch, Disp   = 0; end

%-Time and Frequency windows of interest
%--------------------------------------------------------------------------
woi          = D.inv{val}.contrast.woi;
foi          = D.inv{val}.contrast.fboi;
Nw           = size(woi,1);

%-Get output image field of view and resolution
%--------------------------------------------------------------------------
if space
    sMRIfile = fullfile(spm('dir'),'canonical','avg152T2.nii');
else
    sMRIfile = D.inv{val}.mesh.sMRI;
end
Vin          = spm_vol(sMRIfile);

%-Get surface mesh
%--------------------------------------------------------------------------
if space
    m        = export(gifti(D.inv{val}.mesh.tess_mni),'patch');
else
    m        = export(gifti(D.inv{val}.mesh.tess_ctx),'patch');
end
nd           = D.inv{val}.inverse.Nd;

%-Accumulate mean of log-contrasts (over trials)
%==========================================================================
GL      = spm_mesh_smooth(m);

GW      = D.inv{val}.contrast.GW;

bytrial = iscell(GW{1});

if bytrial
    Ne  = cellfun(@numel,GW);
else
    Ne  = ones(1, numel(GW));
end

Nj      = numel(GW)/Nw;

k  = 1;
iw = [];
ie = [];
for c = 1:length(GW)
    if bytrial
        cGW = GW{c};
    else
        cGW = GW(c);
    end

    for t = 1:Ne(c)
        %-Smooth on the cortical surface
        %------------------------------------------------------------------
        ssq{k} = full(sparse(D.inv{val}.inverse.Is,1,cGW{t},nd,1));
        ssq{k} = spm_mesh_smooth(GL,ssq{k},smooth);

        %-Compute (truncated) moment
        %------------------------------------------------------------------
        lss        = log(ssq{k} + eps);
        i          = lss > (max(lss) - log(32));
        meanlss(k) = mean(lss(i));

        iw(k) = c;
        ie(k) = t;

        k = k + 1;
    end
end

scale = exp(mean(meanlss));

%-Save as meshes
%==========================================================================
if strcmpi(fmt,'mesh')
    [pth,name] = fileparts(D.fname);
    tag = cell(1,Nw);
    for i = 1:Nw
        tag{i} = ['t' sprintf('%g_', woi(i,:)) 'f' sprintf('%d_', foi)];
    end

    %-Save mesh topology
    %----------------------------------------------------------------------
    g = gifti(m);
    save(g,fullfile(D.path,[name '.surf.gii']));

    %-Save mesh data
    %----------------------------------------------------------------------
    spm_progress_bar('Init',numel(ssq),'Exporting as meshes','');
    fprintf('%-40s: %30s','Export as meshes','...please wait');         %-#
    for c = 1:numel(ssq)

        fprintf('%s%30s',repmat(sprintf('\b'),1,30),...
            sprintf('...mesh %d/%d',c,numel(ssq)));                     %-#

        %-Filename
        %------------------------------------------------------------------
        con       = mod(iw(c) - 1, Nj) + 1;
        str       = tag{ceil(iw(c)/Nj)};
        if bytrial, bt = sprintf('_%.0f',ie(c)); else bt = ''; end
        fname     = fullfile(D.path,...
            sprintf('%s_%.0f_%s%.0f%s.gii', name, val, str, con, bt));

        %-Normalise
        %------------------------------------------------------------------
        Contrast = ssq{c} / scale;
        %Contrast = Contrast.*(Contrast > exp(-8));

        %-Write mesh
        %------------------------------------------------------------------
        g = gifti;
        g.cdata = Contrast;
        g.private.metadata(1).name = 'SurfaceID';
        g.private.metadata(1).value = [name '.surf.gii'];
        save(g, fname, 'ExternalFileBinary');

        %-Store filename
        %------------------------------------------------------------------
        D.inv{val}.contrast.fname{c} = fname;

        spm_progress_bar('Set', c);

    end

    spm_progress_bar('Clear');
    fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done');           %-#

    return;
end

%-Normalise and embed in 3D-space
%==========================================================================
fprintf('%-40s: %30s','Writing images','...please wait');               %-#

[pth,name] = fileparts(D.fname);
tag = cell(1,Nw);
for i = 1:Nw
    tag{i} = ['t' sprintf('%d_', woi(i,:)) 'f' sprintf('%d_', foi)];
end

spm_progress_bar('Init',numel(ssq),'Interpolating images','');

con       = mod(iw - 1, Nj) + 1;
win       = ceil(iw/Nj);

ucon      = unique(con);
uwin      = unique(win);

n = 0;
for c = 1:length(ucon)
    for w = 1:numel(uwin)

        ind = find((con == ucon(c)) & (win == uwin(w)));
        str = tag{uwin(w)};

        fname     = fullfile(D.path,...
            sprintf('%s_%.0f_%s%.0f.nii', name, val, str, ucon(c)));

        %-Initialise image
        %------------------------------------------------------------------
        N      = nifti;
        N.dat  = file_array(fname, [Vin.dim, length(ind)], 'FLOAT32-LE');
        N.mat  = Vin.mat;
        N.mat0 = Vin.mat;
        create(N);

        for i = 1:length(ind)

            n = n+1;

            fprintf('%s%30s',repmat(sprintf('\b'),1,30),...
                sprintf('...image %d/%d',n,numel(ssq)));                %-#

            %-Normalise
            %--------------------------------------------------------------
            Contrast = ssq{ind(i)} / scale;

            %-Interpolate those values into voxels
            %--------------------------------------------------------------
            RECimage = spm_mesh_to_grid(m, Vin, Contrast);

            %-3D smoothing and thresholding
            %--------------------------------------------------------------
            spm_smooth(RECimage, RECimage, 1);
            RECimage = RECimage.*(RECimage > exp(-8));

            %-Write (smoothed and scaled) image
            %--------------------------------------------------------------
            N.dat(:, :, :, i) = RECimage;

            %-Store filename
            %--------------------------------------------------------------
            D.inv{val}.contrast.fname{n} = fname;

            spm_progress_bar('Set', n);

        end
    end
end

spm_progress_bar('Clear');
fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done');               %-#

%-Display
%==========================================================================
if Disp && ~spm('CmdLine'), spm_eeg_inv_image_display(D); end

author: spm


text:
function [net, opts, in_name, in_dim] = resnet50(opts)
opts.imageSize = 224;
opts.maxGpuImgs = 90;  % # of images that a 12GB GPU can hold

net = load(fullfile(opts.localDir, 'models', 'imagenet-resnet-50-dag.mat'));
net = dagnn.DagNN.loadobj(net) ;

% remove softmax, fc
net.removeLayer('prob');
net.removeLayer('fc1000');

% freeze pretrained layers
if opts.lastLayer
    for i = 1:numel(net.params)
        net.params(i).learningRate = 0;
        net.params(i).weightDecay  = 0;
    end
end

in_name = 'pool5';
in_dim = 2048;
end

author: kunhe

repo: FastAP-metric-learning


text:
function [shape_initial] = resetshape(bbox, shape_union)
%RESETSHAPE Summary of this function goes here
%   Function: reset the initial shape according to the groundtruth shape and union shape for all faces
%   Detailed explanation goes here
%   Input:
%       bbox: bbounding box of groundtruth shape
%       shape_union: uniionshape
%   Output:
%       shape_initial: reset initial shape
%       bbox: bounding box of face image

% get the bounding box according to the ground truth shape
width_union = (max(shape_union(:, 1)) - min(shape_union(:, 1)));
height_union = (max(shape_union(:, 2)) - min(shape_union(:, 2)));

shape_union = bsxfun(@minus, (shape_union), (min(shape_union)));

shape_initial = bsxfun(@times, shape_union, [(bbox(3)/width_union) (bbox(4)/height_union)]);
shape_initial = bsxfun(@plus, shape_initial, double([bbox(1) bbox(2)]));

end

author: jwyang

repo: face-alignment


text:
function [M,delta_list,theta_list] = compute_edge_patches(w, options)

% [M,delta_list,theta_list] = compute_edge_patches(w, options);

options.null = 0;
if isfield(options, 'n_delta')
    n_delta = options.n_delta;
else
    n_delta = 11;
end
if isfield(options, 'n_theta')
    n_theta = options.n_theta;
else
    n_theta = 12;
end
if isfield(options, 'sigma')
    sigma = options.sigma;
else
    sigma = 0.1;
end
sigma = sigma*w;
if isfield(options, 'rescale')
    rescale = options.rescale;
else
    rescale = 1;
end
% the window to weight the distance
if isfield(options, 'bell')
    bell = options.bell;
else
    bell = 'constant';
end
if isfield(options, 'manifold_type')
    manifold_type = options.manifold_type;
else
    manifold_type = 'edges';
end

switch lower(bell)
    case 'constant'
        B = ones(w);
    case 'sine'
        x = linspace(0,1,w);
        x = (1-cos(2*pi*x))/2;
        B = x'*x;
    otherwise
        error('Unknown bell shape');
end

dmax = sqrt(2)*w/2+2*sigma;
eta = 2;
t = linspace(0,2*pi-2*pi/n_theta,n_theta); t = t(:);
d = linspace(-1,1,n_delta);
d = sign(d).*abs(d).^eta;
d = d(:)*dmax;

[theta_list,delta_list] = meshgrid( t, d );
theta_list = theta_list(:);
delta_list = delta_list(:);
p = length(delta_list);

x = linspace( -(w-1)/2,(w-1)/2, w );
[Y,X] = meshgrid(x,x);
M = zeros(w,w,p);
% compute images
for k=1:p
    t = theta_list(k);
    d = delta_list(k);
    y = cos(t)*X + sin(t)*Y - d;
    if strcmp(manifold_type, 'edges')
        A = tanh(y/sigma);
    else
        A = 2*( 1 - exp( -y.^2 / (2*sigma^2) ) ) - 1;
    end
    M(:,:,k) = A;
end

% normalize
if rescale
    B = repmat(B,[1,1,p]);
    s = sqrt( sum( sum(B .* (M.^2),1), 2) );
    M = M ./ repmat(s,[w,w,1]);
end

author: gpeyre

repo: matlab-toolboxes


text:
function C = mtimes(A,B)
%MTIMES tensor-scalar multiplication.
%
%   C = MTIMES(A,B) is called for the syntax 'A * B' when A or B is a
%   tensor and the other argument is a scalar.
%
%   For tensor-matrix multiplication, use TTM.
%   For tensor-tensor multiplication, use TTT.
%   For tensor-tensor array multiplication, use TIMES or 'A .* B'.
%
%   Examples
%   X = tenrand([3,4,2])
%   W = 5 * X
%
%   See also TENSOR, TENSOR/TTM, TENSOR/TTT, TENSOR/TIMES
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%%
if isscalar(B)
    C = A;
    C.data = B * C.data;
    return;
end

if isscalar(A)
    C = B;
    C.data = A * C.data;
    return;
end

error('Mtimes only supports a tensor times a scalar');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%SOLVE  Solves one or more linear systems or least-squares problems
%
%     [dst,ret] = cv.solve(src1, src2)
%     [...] = cv.solve(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __src1__ input matrix on the left-hand side of the system.
% * __src2__ input matrix on the right-hand side of the system.
%
% ## Output
% * __dst__ output solution.
% * __ret__ Logical return value, see below.
%
% ## Options
% * __Method__ solution (matrix inversion) method, default 'LU'. One of the
%   following matrix decomposition types:
%   * __LU__ Gaussian elimination with the optimal pivot element chosen.
%   * __SVD__ singular value decomposition (SVD) method; the system can be
%     over-defined and/or the matrix `src1` can be singular.
%   * __EIG__ eigenvalue decomposition; the matrix `src1` must be symmetrical.
%   * __Cholesky__ Cholesky LLT factorization; the matrix `src1` must be
%     symmetrical and positively defined.
%   * __QR__ QR factorization; the system can be over-defined and/or the
%     matrix `src1` can be singular.
% * __IsNormal__ this flag can be used together with any of the previous
%   methods; it means that the normal equations
%   `src1' * src1 * dst = src1' * src2` are solved instead of the original
%   system `src1 * dst = src2`. defaul false
%
% The function cv.solve solves a linear system or least-squares problem (the
% latter is possible with 'SVD' or 'QR' methods, or by specifying the flag
% 'IsNormal'):
%
%     dst = argmin_{X} ||src1*X - src2||
%
% If 'LU' or 'Cholesky' method is used, the function returns 1 if `src1`
% (or `src1'*src1`) is non-singular. Otherwise, it returns 0. In the latter
% case, `dst` is not valid. Other methods find a pseudo-solution in case of a
% singular left-hand side part.
%
% ### Note
% If you want to find a unity-norm solution of an under-defined singular
% system `src1 * dst = 0`, the function cv.solve will not do the work.
% Use cv.SVD.SolveZ instead.
%
% See also: cv.invert, cv.SVD, cv.eigen, mldivide, inv, svd, eig
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% VL_DEMO_SIFT_OR  Demonstrates SIFT orientation detection

I = vl_impattern('wedge') ;
ur = 1:size(I,2) ;
vr = 1:size(I,1) ;

% distribute frames on a grid
[u,v] = meshgrid(ur(5:10:end-4),vr(5:10:end-4)) ;
f = [u(:)';v(:)'] ;
K = size(f,2) ;
f = [f ; 4 * ones(1,K) ; 0 * ones(1,K)] ;

% detect orienntations
f = vl_sift(single(I), 'frames', f, 'orientations') ;

figure(1) ; clf ;
imagesc(single(I)) ; colormap gray ; hold on ;
vl_plotframe(f,'color','k','linewidth',3) ;
vl_plotframe(f,'color','y','linewidth',2) ;
axis equal ; axis off ;
vl_demo_print('sift_or') ;

author: jianxiongxiao

repo: ProfXkit


text:
function [zProj,AProj]=projEllipsOntoSpace(z,A,gammaVal,d,T)
%%PROJELLIPSONTOSPACE Determine the orthogonal projection of an ellipsoid
%                    onto an affine space. The projection itself will be an
%                    ellipsoid. For example, in 3D, the projection onto a
%                    plane will be an ellipse.
%
%INPUTS: z  The numDimX1 center of the ellipsoid.
%        A  A numDimXnumDim symmetric, positive definite matrix that
%           specifies the size and shape of the ellipse or ellipsoid, where
%           a point zp is on the ellipse/ellipsoid if
%           (zp-z)'*A*(zp-z)=gammaVal.
%  gammaVal The threshold for declaring a point to be in the ellipsoid. If
%           an empty matrix is passed, the default value of 1 is used.
%      d,T  A numDimX1 point d and a numDImXn matrix such that points in
%           the affine subspace are defined by x=d+T*t where t is an nX1
%           real vector, n<numDim.
%
%OUTPUTS: zProj, AProj Parameters defining the projected ellipsoid. zProj
%           is numDImX1 and AProj is numDimXnumDim. However, AProj will be
%           a singular matrix as everything resides on a subspace. A point
%           is on the projected ellipsoid if
%           (zp-zProj)'*AProj*(zp-zProj)=1
%           The example below shows gives an example of using such a
%           singular matrix to display the projection of a 3D ellipsoid
%           onto a 2D plane.
%
%The algorithm is taken from Section 13 of [1].
%
%EXAMPLE:
%We will project an ellipsoid in 3D onto a plane.
%Define the ellipsoid:
% A=[10,0,0;
%     0,20,0;
%     0,0,1];
% %A rotation matrix
% M=[-0.039622269933159, -0.378261925022885, -0.924850253718582;
%     0.765018603761816, -0.606901366847761,  0.215446668149313;
%    -0.642788154545232, -0.698991163746939,  0.313424219517309];
% A=M*A*M';
% c=[0;0;0];%The center of the ellipsoid
% %Define a plane.
% d=[1;0;3];
% T=[1,1;
%    0,1;
%    0,1];
% T=orth(T);%Must be othonormal.
% %Determine the projection of the ellipsoid onto the plane.
% [cProj,AProj]=projEllipsOntoSpace(c,A,1,d,T);
%
% figure(1)
% clf
% hold on
% axis([-4,4,-4,4,-4,4])
% %Draw the ellipsoid
% drawEllipse(c,A,1,'EdgeColor','none')
% view(104,-29)
% %Compute points on the plane to display.
% numPoints=50;
% tSpan=linspace(-10,10,numPoints);
% [tx,ty]=meshgrid(tSpan,tSpan);
% t=[tx(:)';ty(:)'];
% placeSurfVals=bsxfun(@plus,d,T*t);
% X=reshape(placeSurfVals(1,:),numPoints,numPoints);
% Y=reshape(placeSurfVals(2,:),numPoints,numPoints);
% Z=reshape(placeSurfVals(3,:),numPoints,numPoints);
% %Draw the plane
% surf(X,Y,Z)
% %Draw the projection of the ellipsoid onto the plane
% drawEllipse(cProj,AProj,1)
% %This makes the ellipsoid shiny
% light('Position',[0;1;-1])
%
%REFERENCES:
%[1] S. B. Pope, "Algorithms for ellipsoids," Cornell University, Tech.
%    Rep. FDA-08-01, Feb. 2008. [Online].
%    Available: https://tcg.mae.cornell.edu/pubs/Pope_FDA_08.pdf
%
%May 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(isempty(gammaVal))
   gammaVal=1;
end

A=A/gammaVal;

L=chol(A,'lower');

%Make T into an orthonormal matrix. This does not rotate or shift the
%affine space at all. It just changes the parameterization as in x=d+T*t
%still holds, but if T is orthonormalized, then a different value of t
%is needed to get the same x.
T=orth(T);%Must be othonormal.

%Extracted from Equation 126
zProj=d+T*T'*(z-d);
[U,Sigma,~]=svd(T'/(L'),'econ');

L=T*U/Sigma;
AProj=L*L';

end

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%% Neural Art with lbfgs optimization
%% It cannot run on my laptop T^T. The error seems to be memory out.
caffe.reset_all();
caffe.set_mode_gpu();
gpu_id = 1;  % we will use the first gpu in this demo
caffe.set_device(gpu_id);

% original_file = 'D:\deepLearning\caffe-windows\matlab\demo\Inceptionism\googlenet_neuralart.prototxt';
% net_weights = 'D:\deeplearning\caffe-windows\examples\GoogLeNet\imagenet_googlenet.caffemodel';
% style_layer = {'icp2_in','icp3_out','icp5_out','icp7_out','icp9_out'};
original_file = 'D:\deepLearning\caffe-windows\matlab\demo\NeuralArt\VGG_16_nueralart.prototxt';
net_weights = 'D:\deepLearning\caffe-windows\matlab\demo\NeuralArt\VGG16_thinned_net.caffemodel';
style_layer = {'conv1_1','conv2_1','conv3_1','conv4_1','conv5_1'};
content_layer = {'conv4_2'};

style_weights = [1 1 1 1 1 0.05];
tv_norm = 0.001;
use_color_prior = false;
num_cluster=6;
color_prior = 0.5;
long_size = 512;

style_image = imread('d:\starry_night.jpg');
if size(style_image,1) > size(style_image,2)
    style_image = imresize(style_image,[long_size, size(style_image,2) / size(style_image,1) * long_size]);
else
    style_image = imresize(style_image,[size(style_image,1) / size(style_image,2) * long_size, long_size]);
end;
content_image = imread('d:\hoovertowernight.jpg');
if size(content_image,1) > size(content_image,2)
    content_image = imresize(content_image,[long_size, size(content_image,2) / size(content_image,1) * long_size]);
else
    content_image = imresize(content_image,[size(content_image,1) / size(content_image,2) * long_size, long_size]);
end;
figure(1);
imshow(style_image);
title('style image');
figure(2);
imshow(content_image);
title('content image');
[style_generate_prototxt, style_pattern, content_pattern,colorObj] = MakeStylePrototxt(original_file, net_weights, style_layer, style_weights, content_layer, style_image, content_image,num_cluster);
if use_color_prior
colorObj
num_cluster = colorObj.NComponents;
end;
if use_color_prior
    gaussian_net = caffe.Net('gaussian_net.prototxt','test');

    W = single(zeros(1,1,3,3*num_cluster));
    b = single(zeros(3*num_cluster,1));

    for i=1:num_cluster
        W(:,:,:,(i-1)*3+1:i*3) = inv(colorObj.Sigma(:,:,i));
        b((i-1)*3+1:i*3) = -1 * colorObj.Sigma(:,:,i) \ colorObj.mu(i,:)';
    end;

    nth_layer = gaussian_net.layer_vec(gaussian_net.name2layer_index('gaussian_prior'));
    nth_layer.params(1).set_data(W);
    nth_layer.params(2).set_data(b);
end;

forward_input = cell(length(style_pattern)+2,1);
forward_input(3:end) = style_pattern;
forward_input(2) = {content_pattern};

mean_file = [];
vgg_mean =  [103.939, 116.779, 123.68];
[height, width, channel] = size(content_image);

%%%%%%%%%extract the train features
stylegen_net = caffe.Net(style_generate_prototxt,net_weights,'test');
mean_image = permute(repmat(vgg_mean',[1,size(content_image,2),size(content_image,1)]),[2,3,1]);
input_data = randn(size(mean_image,1), size(mean_image,2), 3, 1, 'single')*50;

lr = 10;
max_lr = 50;
momentum = 0.8;
lastgrad = zeros(size(mean_image));
last_cost = 9999999999999;

for iter = 1:2000
    bak_data = input_data;
    bak_grad = lastgrad;
    [cost, grad] = NeuralArtCost(input_data,[],stylegen_net,forward_input,style_weights);

    if tv_norm > 0
        I = input_data(:,:,:,1);
        Gx = (I(2:end-1,:,:) - I(1:end-2,:,:)) - (I(3:end,:,:) - I(2:end-1,:,:));
        Gx = [(I(1,:,:) - I(2,:,:)); Gx; (I(end,:,:) - I(end-1,:,:))];
        Gy = (I(:,2:end-1,:) - I(:,1:end-2,:)) - (I(:,3:end,:) - I(:,2:end-1,:));
        Gy = [(I(:,1,:) - I(:,2,:)) Gy (I(:,end,:) - I(:,end-1,:))];
        grad = grad + tv_norm * (Gx+Gy);
    end;

    if use_color_prior
        gmm_prior = gaussian_net.forward({input_data});
        sum_gp = zeros(size(mean_image,1),size(mean_image,2));
        sum_prob_gradient = zeros(size(mean_image));
        for i=1:num_cluster
            gp = bsxfun(@minus,input_data(:,:,:,1),reshape(colorObj.mu(i,:),[1 1 3])) .* gmm_prior{1}(:,:,(i-1)*3+1:i*3);
            gp = sum(gp,3);
            gp = colorObj.PComponents(i) * exp(-gp);
            sum_prob_gradient = sum_prob_gradient + bsxfun(@times,gp,gmm_prior{1}(:,:,(i-1)*3+1:i*3));
            sum_gp = sum_gp + gp;
        end;
        sum_prob_gradient = bsxfun(@rdivide,sum_prob_gradient,sum_gp);
        sum_prob_gradient(isnan(sum_prob_gradient)) = 0;
        input_data(:,:,:,1) = input_data(:,:,:,1) - lr * color_prior * sum_prob_gradient;
    end;
    %%%%%%%%%%%%%%%%%%%%%%%%gd linear search
    lastgrad = (1 - momentum) * lr * grad   + momentum * lastgrad;%/ norm(res(:))
    input_data(:,:,:,1) = input_data(:,:,:,1) - lastgrad;
    if cost>last_cost
        lr = lr * 0.9;
%         input_data = bak_data;%why...
%         last_grad = bak_grad;
    end;
    if cost<last_cost&&lr<max_lr%&& (this_prob-last_prob) / this_prob < 0.001
        lr = lr * 1.1;
    end;
    if cost<last_cost
        last_cost = cost;
    end;
    fprintf('iter=%d,lr=%f,this_cost=%f,last_cost=%f\n',iter,lr,cost, last_cost);
    %%%%%%%%%%%%%%%%%%%%%%gd

    if mod(iter,10)==0
        figure(3);
        output = mean_image + input_data(:,:,:,1);
        output = output(:, :, [3, 2, 1]);
        output = permute(output, [2 1 3]);
        imshow(uint8(output));
        I = output;
        title('generated image');
    end;
end;

author: happynear

repo: DeepVisualization


text:
function plot_accuracy_comparison

color = {'r', 'g'};
leng = {'PoseCNN', 'PoseCNN New'};
aps = zeros(2, 1);
lengs = cell(2, 1);
close all;

% load results
object = load('results_comparison.mat');
distances_sys = object.distances_sys;
distances_non = object.distances_non;
rotations = object.errors_rotation;
translations = object.errors_translation;
cls_ids = object.results_cls_id;

index_plot = [1, 2];

% read class names
fid = fopen('classes.txt', 'r');
C = textscan(fid, '%s');
classes = C{1};
classes{end+1} = 'All 21 objects';
fclose(fid);

hf = figure('units','normalized','outerposition',[0 0 1 1]);
font_size = 24;
max_distance = 0.1;

% for each class
for k = 1:numel(classes)
    index = find(cls_ids == k);
    if isempty(index)
        index = 1:size(distances_sys,1);
    end

    % distance symmetry
    subplot(2, 2, 1);
    for i = index_plot
        D = distances_sys(index, i);
        D(D > max_distance) = inf;
        d = sort(D);
        n = numel(d);
        accuracy = cumsum(ones(1, n)) / n;
        plot(d, accuracy, color{i}, 'LineWidth', 4);
        aps(i) = VOCap(d, accuracy);
        lengs{i} = sprintf('%s (%.2f)', leng{i}, aps(i) * 100);
        hold on;
    end
    hold off;
    %h = legend('network', 'refine tranlation only', 'icp', 'stereo translation only', 'stereo full', '3d coordinate');
    %set(h, 'FontSize', 16);
    h = legend(lengs(index_plot), 'Location', 'southeast');
    set(h, 'FontSize', font_size);
    h = xlabel('Average distance threshold in meter (symmetry)');
    set(h, 'FontSize', font_size);
    h = ylabel('accuracy');
    set(h, 'FontSize', font_size);
    h = title(classes{k}, 'Interpreter', 'none');
    set(h, 'FontSize', font_size);
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', font_size)

    % distance non-symmetry
    subplot(2, 2, 2);
    for i = index_plot
        D = distances_non(index, i);
        D(D > max_distance) = inf;
        d = sort(D);
        n = numel(d);
        accuracy = cumsum(ones(1, n)) / n;
        plot(d, accuracy, color{i}, 'LineWidth', 4);
        aps(i) = VOCap(d, accuracy);
        lengs{i} = sprintf('%s (%.2f)', leng{i}, aps(i) * 100);
        hold on;
    end
    hold off;
    %h = legend('network', 'refine tranlation only', 'icp', 'stereo translation only', 'stereo full', '3d coordinate');
    %set(h, 'FontSize', 16);
    h = legend(lengs(index_plot), 'Location', 'southeast');
    set(h, 'FontSize', font_size);
    h = xlabel('Average distance threshold in meter (non-symmetry)');
    set(h, 'FontSize', font_size);
    h = ylabel('accuracy');
    set(h, 'FontSize', font_size);
    h = title(classes{k}, 'Interpreter', 'none');
    set(h, 'FontSize', font_size);
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', font_size)

    % rotation
    subplot(2, 2, 3);
    for i = index_plot
        D = rotations(index, i);
        d = sort(D);
        n = numel(d);
        accuracy = cumsum(ones(1, n)) / n;
        plot(d, accuracy, color{i}, 'LineWidth', 4);
        hold on;
    end
    hold off;
    %h = legend('network', 'refine tranlation only', 'icp', 'stereo translation only', 'stereo full', '3d coordinate');
    %set(h, 'FontSize', 16);
    h = legend(leng(index_plot), 'Location', 'southeast');
    set(h, 'FontSize', font_size);
    h = xlabel('Rotation angle threshold');
    set(h, 'FontSize', font_size);
    h = ylabel('accuracy');
    set(h, 'FontSize', font_size);
    h = title(classes{k}, 'Interpreter', 'none');
    set(h, 'FontSize', font_size);
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', font_size)

    % translation
    subplot(2, 2, 4);
    for i = index_plot
        D = translations(index, i);
        D(D > max_distance) = inf;
        d = sort(D);
        n = numel(d);
        accuracy = cumsum(ones(1, n)) / n;
        plot(d, accuracy, color{i}, 'LineWidth', 4);
        hold on;
    end
    hold off;
    h = legend(leng(index_plot), 'Location', 'southeast');
    set(h, 'FontSize', font_size);
    h = xlabel('Translation threshold in meter');
    set(h, 'FontSize', font_size);
    h = ylabel('accuracy');
    set(h, 'FontSize', font_size);
    h = title(classes{k}, 'Interpreter', 'none');
    set(h, 'FontSize', font_size);
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', font_size)

    filename = sprintf('plots/%s.png', classes{k});
    hgexport(hf, filename, hgexport('factorystyle'), 'Format', 'png');
end

function ap = VOCap(rec, prec)

index = isfinite(rec);
rec = rec(index);
prec = prec(index)';

mrec=[0 ; rec ; 0.1];
mpre=[0 ; prec ; prec(end)];
for i = 2:numel(mpre)
    mpre(i) = max(mpre(i), mpre(i-1));
end
i = find(mrec(2:end) ~= mrec(1:end-1)) + 1;
ap = sum((mrec(i) - mrec(i-1)) .* mpre(i)) * 10;

author: yuxng


text:
%% Canny Edge Detection
% This sample demonstrates Canny edge detection.
%
% In this demo, we show how to use the OpenCV function |cv.Canny| to implement
% the Canny Edge Detector.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/edge.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/edge.py>
% * <https://docs.opencv.org/3.2.0/da/d5c/tutorial_canny_detector.html>
% * <https://docs.opencv.org/3.2.0/da/d22/tutorial_py_canny.html>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/ImgTrans/CannyDetector_Demo.cpp>
%

%% Theory
%
% The _Canny Edge detector_ was developed by John F. Canny in 1986. Also known
% to many as the _optimal detector_, the Canny algorithm aims to satisfy three
% main criteria:
%
% * *Low error rate:* Meaning a good detection of only existent edges.
% * *Good localization:* The distance between edge pixels detected and real
%   edge pixels have to be minimized.
% * *Minimal response:* Only one detector response per edge.
%
%% Steps
%
% 1) Filter out any noise. The Gaussian filter is used for this purpose. An
%    example of a Gaussian kernel of $size = 5$ that might be used is shown
%    below:
%
% $$K = \frac{1}{159}
%       \left[{\matrix{
%           2 &  4 &  5 &  4 & 2 \cr
%           4 &  9 & 12 &  9 & 4 \cr
%           5 & 12 & 15 & 12 & 5 \cr
%           4 &  9 & 12 &  9 & 4 \cr
%           2 &  4 &  5 &  4 & 2
%       }}\right]$$
%
% 2) Find the intensity gradient of the image. For this, we follow a procedure
%    analogous to Sobel:
%
% * Apply a pair of convolution masks (in $x$ and $y$ directions):
%
% $$G_{x} = \left[{\matrix{
%               -1 & 0 & +1 \cr
%               -2 & 0 & +2 \cr
%               -1 & 0 & +1
%           }}\right]$$
%
% $$G_{y} = \left[{\matrix{
%               -1 & -2 & -1 \cr
%                0 &  0 &  0 \cr
%               +1 & +2 & +1
%           }}\right]$$
%
% * Find the gradient strength and direction. The direction is rounded to
%   one of four possible angles (namely 0, 45, 90 or 135)
%
% $$G = \sqrt{ G_{x}^{2} + G_{y}^{2} }$$
%
% $$\theta = \arctan(\frac{ G_{y} }{ G_{x} })$$
%
% 3) _Non-maximum suppression_ is applied. This removes pixels that are not
%    considered to be part of an edge. Hence, only thin lines (candidate edges)
%    will remain.
%
% 4) _Hysteresis_: The final step. Canny does use two thresholds
%    (_upper_ and _lower_). Canny recommended a |upper:lower| ratio between
%    2:1 and 3:1.
%
% * If a pixel gradient is higher than the _upper_ threshold, the pixel is
%   accepted as an edge
% * If a pixel gradient value is below the _lower_ threshold, then it is
%   rejected.
% * If the pixel gradient is between the two thresholds, then it will be
%   accepted only if it is connected to a pixel that is above the _upper_
%   threshold.
%
% For more details, you can always consult your favorite Computer Vision book.
%

%% Code
%
% This program:
%
% * Asks the user to enter a numerical value to set the lower threshold for
%   our _Canny Edge Detector_ (by means of a slider)
% * Applies the _Canny Detector_ and generates a *mask* (bright lines
%   representing the edges on a black background)
% * Applies the mask obtained on the original image and display it in a window
%

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

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

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

    % retrieve current values from UI controls
    apertures = [3, 5, 7];
    aIdx = get(h.pop, 'Value');
    thresh = round(get(h.slid, 'Value'));
    set(h.txt, 'String',sprintf('Threshold: %3d',thresh));

    % convert image to grayscale, and blur to reduce the noise
    gray = cv.cvtColor(h.src, 'RGB2GRAY');
    gray = cv.blur(gray, 'KSize',[3 3]);

    % detect edges, with 3:1 as threshold ratio
    if true
        % default canny (Sobel gradient)
        edges = cv.Canny(gray, thresh*[1 3], 'ApertureSize',apertures(aIdx));
    else
        % canny with custom gradient (Scharr)
        dx = cv.Scharr(gray, 'DDepth','int16', 'XOrder',1, 'YOrder',0);
        dy = cv.Scharr(gray, 'DDepth','int16', 'XOrder',0, 'YOrder',1);
        edges = cv.Canny2(dx, dy, thresh*[1 3]);
    end

    % apply edges mask on original image
    if true
        out = cv.copyTo(h.src, 'Mask',edges);
    else
        out = bsxfun(@times, h.src, uint8(edges~=0));
    end

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

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

    % parameters
    thresh = 10;
    max_thresh = 150;
    sz = size(img);
    sz(2) = max(sz(2), 300);  % minimum figure width

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

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


text:
function y = vl_nnconcat(inputs, dim, dzdy, varargin)
%VL_NNCONCAT CNN concatenate multiple inputs.
%  Y = VL_NNCONCAT(INPUTS, DIM) concatenates the inputs in the cell
%  array INPUTS along dimension DIM generating an output Y.
%
%  DZDINPUTS = VL_NNCONCAT(INPUTS, DIM, DZDY) computes the derivatives
%  of the block projected onto DZDY. DZDINPUTS has one element for
%  each element of INPUTS, each of which is an array that has the same
%  dimensions of the corresponding array in INPUTS.

% Copyright (C) 2015 Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.inputSizes = [] ;
opts = vl_argparse(opts, varargin, 'nonrecursive') ;

if nargin < 2, dim = 3; end;
if nargin < 3, dzdy = []; end;

if isempty(dzdy)
  y = cat(dim, inputs{:});
else
  if isempty(opts.inputSizes)
    opts.inputSizes = cellfun(@size, inputs, 'UniformOutput', false) ;
  end
  start = 1 ;
  y = cell(1, numel(opts.inputSizes)) ;
  s.type = '()' ;
  s.subs = {':', ':', ':', ':'} ;
  for i = 1:numel(opts.inputSizes)
    stop = start + opts.inputSizes{i}(dim) ;
    s.subs{dim} = start:stop-1 ; ;
    y{i} = subsref(dzdy,s) ;
    start = stop ;
  end
end

repo: CREST-Release


text:
function [ y2, m2, d2, f2 ] = yjf_to_ymdf_julian ( y1, j1, f1 )

%*****************************************************************************80
%
%% YJF_TO_YMDF_JULIAN converts a YJF to YMDF date, both in the Julian calendar.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%   13 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, J1, real F1, the YJF date.
%
%    Output, integer Y2, M2, D2, real F2,
%    the YMDF date.
%

%
%  Copy the input.
%
  y2 = y1;
  j2 = j1;
  f2 = f1;
%
%  Check the input.
%
  [ y2, j2, ierror ] = yj_check_julian ( y2, j2 );

  if ( ierror ~= 0 )
    y2 = 0;
    m2 = 0;
    d2 = 0;
    f2 = 0.0;
    return
  end
%
%  Convert the input.
%
  m2 = 1;
  d2 = j2;

  [ y2, m2, d2 ] = day_borrow_julian ( y2, m2, d2 );

  [ y2, m2, d2 ] = day_carry_julian ( y2, m2, d2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function cout = medianfilt(cin, len)

%tstoolbox/@core/medianfilt
%   Syntax:
%     * medianfilt(cin,len)
%
%   Input Arguments:
%     * cin - core object
%
%   moving median filter
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

N = dlens(cin,1); 	% Laenge in der ersten Dimension (sollte auch die einzige sein (skalar))
ts = data(cin);		% Die Datenwerte des Input-Signals werden nach ts kopiert

md = zeros(N-len+1,1);

for i=1:N-len+1
	md(i) = median(ts(i:i+len-1));
end

cout = core(md);