text:
function [k, rbfPart, n2] = sqexpKernCompute(kern, x, x2)

% SQEXPKERNCOMPUTE Compute the SQEXP kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the pre-built compound squared exponential
% kernel given inputs associated with rows and columns.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 : the input matrix associated with the columns of the kernel.
% RETURN k : the kernel matrix computed at the given points.
%
% FORMAT
% DESC computes the kernel matrix for the pre-built compound squared exponential
% kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : the kernel matrix computed at the given points.
%
% SEEALSO : sqexpKernParamInit, kernCompute, kernCreate, sqexpKernDiagCompute
%
% COPYRIGHT : Neil D. Lawrence, 2004

% KERN

if nargin < 3
  n2 = dist2(x, x);
  wi2 = (.5 .* kern.inverseWidth);
  rbfPart = kern.rbfVariance*exp(-n2*wi2);
  k = rbfPart + kern.whiteVariance*eye(size(x, 1));
else
  n2 = dist2(x, x2);
  wi2 = (.5 .* kern.inverseWidth);
  rbfPart = kern.rbfVariance*exp(-n2*wi2);
  k = rbfPart;
end
k = k + kern.biasVariance;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function res=mocov_find_files(root_dir, file_pat, monitor, exclude_pat)
% Finds files recursively in a directory
%
% res=mocov_find_files([root_dir[, file_pat]])
%
% Inputs:
%   root_dir            Optional directory in which files are sought.
%                       If not provided or empty, the current working
%                       directory is used.
%   file_pat            Optional wildcard pattern, e.g. '*.m' for all
%                       files ending with '.m'. If omitted, the pattern '*'
%                       is used, corresponding to all files.
%   monitor             Optional progress monitory that supports a
%                       'notify' method.
%   exclude_pat         Optional cell array of patterns to exclude. Both
%                       files and directories which match one of these
%                       patterns will be omitted from the output.
%
% Output:
%   res                 Kx1 cell with names of files in root_dir matching
%                       the pattern.
%
% NNO May 2015

    if nargin<1
        root_dir='';
    end

    if ~(isempty(root_dir) || mocov_util_isfolder(root_dir))
        error('first argument must be directory');
    end

    if nargin<2
        file_pat='*';
    end

    if nargin<3
        monitor=[];
    end

    if nargin<4
        exclude_pat={};
    end

    if ischar(exclude_pat)
        exclude_pat={exclude_pat};
    end

    file_re=pattern2re(file_pat);
    exclude_re=get_exclude_re(exclude_pat);

    if ~isempty(monitor)
        msg=sprintf('Finding files matching %s from %s',file_pat,root_dir);
        notify(monitor, msg);
    end

    res=find_files_recursively(root_dir,file_re,monitor,exclude_re);

function re=pattern2re(pat)
    re=['^' ... % start of the string
        regexptranslate('wildcard',pat) ...
        '$'];   % end of the string

function re=get_exclude_re(exclude_pat)
    n=numel(exclude_pat);
    if n==0
        re='';
        return;
    end

    excl_re_cell=cellfun(@pattern2re,exclude_pat,...
                        'UniformOutput',false);

    joined=sprintf('|%s',excl_re_cell{:});
    re=joined(2:end);

function res=find_files_recursively(root_dir,file_re,monitor,exclude_re)
    if isempty(root_dir)
        dir_arg={};
    else
        dir_arg={root_dir};
    end

    d=dir(dir_arg{:});
    n=numel(d);

    res_cell=cell(n,1);
    for k=1:n
        fn=d(k).name;

        ignore_fn=strcmp(fn,'.') || strcmp(fn,'..');

        res=cell(0,1);
        if ~ignore_fn
            path_fn=fullfile(root_dir, fn);

            if ~isempty(regexp(fn,exclude_re,'once'));
                continue;
            elseif mocov_util_isfolder(path_fn)
                res=find_files_recursively(path_fn,file_re,...
                                                monitor,exclude_re);
            elseif ~isempty(regexp(fn,file_re,'once'));
                res={path_fn};
                if ~isempty(monitor)
                    notify(monitor,'.',path_fn);
                end
            end
        end
        res_cell{k}=res;
    end

    res_cell=res_cell(~cellfun(@isempty,res_cell));
    res=cat(1,res_cell{:});

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%% 2D Riemann Problem: Dimensional Splitting with Front Tracking
% In this example, we consider the 2D scalar conservation law
%
% $$ u_t + (u^2)_x/2 + (u^2)_y/2 = 0$$
%
% with Riemann initial data, that is, with an intial function that has a
% constant value in each of the four quadrants and discontinuities along
% the x- and y-axes.
%
% To solve the one-dimensional problems arising as part of the operator
% splitting, we will use front tracking followed by a projection step onto
% underlying regular grid.

%% Initial setup
N = 100;
h = 2/N;
x = -1:h:1;
y = 0.5*(x(1:end-1)+x(2:end));
T = 1.0;
delta = 0.1/sqrt(N);
[X,Y]=meshgrid(y,y);

%% Exact solutions
% We consider two different setups. In both setups u0=1/4 in the second and
% fourth quadrant. Setting u0=1/2 in the first quadrant and u0=-1/2 in the
% third quadrant gives a self-similar wave pattern where the four original
% constant states are separated by two pairs of rarefaction waves. For each
% pair, the two rarefaction waves meet in a sharp kink along the line y =
% x. Switching the values in the first and third quadrant gives a wave
% pattern that consists of the four constant states separated by six shocks
% forming two triple points at (3t/8,−t/8) and (−t/8, 3t/8).
subplot(2,2,1), pcolor(X,Y,rpsola(X,Y,0)), shading flat; colorbar, title('Problem A')
subplot(2,2,2), pcolor(X,Y,rpsolb(X,Y,0)), shading flat; colorbar, title('Problem B')
subplot(2,2,3), pcolor(X,Y,rpsola(X,Y,1)), shading flat; colorbar
subplot(2,2,4), pcolor(X,Y,rpsolb(X,Y,1)), shading flat; colorbar

%% Riemann problem B: shocks
% First, we study the error for the second Riemann problem as a function of
% the number of operator splitting steps
u0     = rpsolb(X,Y,0);
utrue  = rpsolb(X,Y,1);
rel    = sum(sum(abs(utrue)));
for i=1:4,
  Nstep = 4^(i-1);
  u     = DimSplit(u0,x,x,'burgerRsol','burgerRsol',delta,Nstep,T,'wbar');
  feil  = sum(sum(abs(utrue-u)))/rel;

  subplot(2,2,i);
  pcolor(y,y,u), axis equal image;  colorbar, shading flat
  title(sprintf('Grid: %dx%d  Steps: %d', N, N, Nstep));
  xlabel(sprintf('L1 error: %4.2f %%', 100*feil));
  p=get(gca,'position'); p([3 4])=p([3 4])+0.03;
  set(gca,'position',p,'XTickLabel',[]);
end;
%%
% From the figure, we see that the error is a convex function of the number
% of splitting steps. With few splitting steps, the error is dominated by
% the splitting error and the strong self-sharpening effects present near
% the shocks will counteract the smearing introduced by the projection
% operator that brings the front tracking solution back from a possibly
% irregular grid and back to the underlying Cartesian grid. As the number
% of time steps increases, the number of projections increases and the
% self-sharpening effects of diminish so that the smearing errors become
% dominant. Because the two error mechanisms work in different directions,
% a minimum error is observed for an intermediate number of splitting
% steps.

%% Convergence test
% We perform a grid-refinement study for the two Riemann problems, fixing
% the CFL number to 10 and report the error as a logarithmic plot
nu=10;
errA=ones(1,6); errB=errA;
for i=1:6,
  N     = 8*2^i;
  h     = 2/N;
  x     = -1:h:1;
  y     = 0.5*(x(1:end-1)+x(2:end));
  delta = 0.1/sqrt(N);
  [X,Y] = meshgrid(y,y);

  % Riemann problem A
  u       = DimSplit(rpsola(X,Y,0), x,x,'burgerRsol','burgerRsol', ...
                     delta, round(T/(nu*h)), T,'wbar');
  utrue   = rpsola(X,Y,1);
  rel     = sum(sum(abs(utrue)));
  errA(i) = sum(sum(abs(utrue-u)))/rel;

  % Riemann problem B
  u       = DimSplit(rpsolb(X,Y,0), x,x,'burgerRsol','burgerRsol', ...
                     delta, round(T/(nu*h)), T,'wbar');
  utrue   = rpsolb(X,Y,1);
  rel     = sum(sum(abs(utrue)));
  errB(i) = sum(sum(abs(utrue-u)))/rel;
end
clf, semilogy(4:9,errA,'-o',4:9,errB,'x-')
xlabel('log2 of grid size'); ylabel('Relative L1 error');
legend('Problem A','Problem B');
%%
% As pointed out above, the strong self-sharpening effects present near the
% shocks in Riemann problem B will counteract the smearing introduced by
% the projection operator leading to a convergence rate of approximately
% one. For the smooth solution resulting from Riemann problem A, the
% self-sharpening effects are much weaker, thereby giving a lower
% convergence rate.

repo: Aero-matlab


text:
function [ rexp, sexp ] = poly_q4 ( )

%*****************************************************************************80
%
%% POLY_Q4 returns the monomials associated with a 4 node quadrilateral.
%
%  Reference Element Q4:
%
%    |
%    1  4-----3
%    |  |     |
%    |  |     |
%    S  |     |
%    |  |     |
%    |  |     |
%    0  1-----2
%    |
%    +--0--R--1-->
%
%  Formula:
%
%    Given coefficients A(I), the polynomial interpolant at (R,S) is
%
%      P(R,S) = sum ( 1 <= I <= N ) A(I) * R**REXP(I) * S**SEXP(I)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer REXP(4), SEXP(4), the powers of R and S associated
%    with each polynomial.
%
  rexp(1:4) = [ 0, 0, 1, 1 ];
  sexp(1:4) = [ 0, 1, 0, 1 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function D = addcroplabel(D, words)
%
% Adds the field object.crop (it also considers an occluded object as a
% cropped object).
%
% Removes the words 'crop', 'occluded', etc. from the object name.

if nargin < 2
    words = {'crop', 'occluded', 'part'};
end

Nwords = length(words);

for i = 1:length(D)
    if isfield(D(i).annotation, 'object')
        Nobjects = length(D(i).annotation.object);
        for j = 1:Nobjects
            if isfield(D(i).annotation.object(j), 'occluded')
                if strcmp(D(i).annotation.object(j).occluded,'yes')
                    D(i).annotation.object(j).crop = '1';
                end
            end

            w = lower(getwords(D(i).annotation.object(j).name));
            k = ismember(w, words);

            if sum(k)>0
                D(i).annotation.object(j).crop = '1';
                % remove words
                for n = 1:Nwords
                    D(i).annotation.object(j).name = strrep(D(i).annotation.object(j).name, words{n}, '');
                end
            else
                D(i).annotation.object(j).crop = '0';
            end
        end
    end
end

% Make sure there are no new empty object names
D = LMvalidobjects(D);

author: CSAILVision

repo: LabelMeToolbox


text:
function tauxovertau = tau_prob ( fh, fhx, ch, chx )

%*****************************************************************************80
%
%% TAUPROB evaluates the derivative of TAU with respect to (X/TAU).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 November 2005
%
%  Parameters:
%
  global alpha1
  global alpha2
  global beta1
  global beta2
  global gamma1
  global gamma2

  cprod = ( ch + alpha1 ) * ( ch + alpha2 );
  fprod = ( fh + beta1 ) * ( fh + beta2 );
  dela = alpha2 - alpha1;
  delb = beta2 - beta1;
  tauxovertau = gamma1 * dela * chx / cprod + gamma2 * fhx * delb / fprod;

  return
end


text:
function [whichScans,ok] = er_selectScans(vw, title)
% [whichScans, ok] = er_selectScans(vw, title);
%
% Like selectScans, provides a prompt to select a subset
% of scans for the current view/data Type, returning an
% array of the scan numbers selected.
%
% I like this a litte better than select scans, b/c it
% gives the annotation for each scan, and makes multi select
% easier.
%
%
% 04/05/04 ras.
global dataTYPES;

if ~exist('title','var')    title = 'Choose scans:';    end

whichSeries = vw.curDataType;
nscans = viewGet(vw, 'numScans');
scanList    = cell(1, nscans);
for i = 1:nscans
    scanList{i} = dataTYPES(whichSeries).scanParams(i).annotation;
    if isempty(scanList{i}), scanList{i} = sprintf('Scan%i',i); end
end

[whichScans, ok] = listdlg('PromptString',title,...
    'ListSize',[300 400],...
    'ListString',scanList,'InitialValue',1 ,'OKString','OK');

% Maybe we should check the return and announce of there is a problem?

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_bug3218

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_rejectartifact ft_rejectvisual

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug3218.mat'), 'data');

% When rejecting artifacts using 'partial', ft_rejectvisual works fine
cfg = [];
cfg.artfctdef.visual.artifact = data.sampleinfo(1, :) + [300 -300];
cfg.artfctdef.reject = 'partial';
dataReject = ft_rejectartifact(cfg, data);

ft_rejectvisual([], dataReject);

% When rejecting artifacts using 'nan', ft_rejectvisual recalculates the
% metric after each trial that is excluded. With this small data set it's
% not really a problem, but it becomes cumbersome with many trials and
% channels.

cfg = [];
cfg.artfctdef.visual.artifact = data.sampleinfo(1, :) + [300 -300];
cfg.artfctdef.reject = 'nan';
dataReject = ft_rejectartifact(cfg, data);

ft_rejectvisual([], dataReject);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
Network ResNeXt50 {
	Layer CONV1 {
		Type: CONV
		Stride { X: 2, Y: 2 }
		Dimensions { K: 64, C: 3, R: 1, S: 1, Y:224, X:224 }
	}

	Layer CONV2_1_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 128, C: 64, R: 1, S: 1, Y: 56, X: 56 }
	}

	Layer CONV2_1_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 4, C: 128, R: 3, S: 3, Y: 56, X: 56 }
	}

	Layer CONV2_1_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 128, R: 1, S: 1, Y: 56, X: 56 }
	}

	Layer CONV2_2_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 128, C: 256, R: 1, S: 1, Y: 56, X: 56 }
	}

	Layer CONV2_2_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 4, C: 128, R: 3, S: 3, Y: 56, X: 56 }
	}

	Layer CONV2_2_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 128, R: 1, S: 1, Y: 56, X: 56 }
	}

	Layer CONV2_3_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 128, C: 256, R: 1, S: 1, Y: 56, X: 56 }
	}

	Layer CONV2_3_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 4, C: 128, R: 3, S: 3, Y: 56, X: 56 }
	}

	Layer CONV2_3_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 128, R: 1, S: 1, Y: 56, X: 56 }
	}

	Layer CONV3_1_1 {
		Type: CONV
		Stride { X: 2, Y: 2 }
		Dimensions { K: 256, C: 256, R: 1, S: 1, Y: 56, X: 56 }
	}

	Layer CONV3_1_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 8, C: 256, R: 3, S: 3, Y: 28, X: 28 }
	}

	Layer CONV3_1_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 256, R: 1, S: 1, Y: 28, X: 28 }
	}

	Layer CONV3_2_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 256, R: 1, S: 1, Y: 28, X: 28 }
	}

	Layer CONV3_2_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 8, C: 256, R: 3, S: 3, Y: 28, X: 28 }
	}

	Layer CONV3_2_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 256, R: 1, S: 1, Y: 28, X: 28 }
	}

	Layer CONV3_3_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 256, R: 1, S: 1, Y: 28, X: 28 }
	}

	Layer CONV3_3_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 8, C: 256, R: 3, S: 3, Y: 28, X: 28 }
	}

	Layer CONV3_3_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 256, R: 1, S: 1, Y: 28, X: 28 }
	}

	Layer CONV3_4_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 256, C: 256, R: 1, S: 1, Y: 28, X: 28 }
	}

	Layer CONV3_4_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 8, C: 256, R: 3, S: 3, Y: 28, X: 28 }
	}

	Layer CONV3_4_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 256, R: 1, S: 1, Y: 28, X: 28 }
	}

	Layer CONV4_1_1 {
		Type: CONV
		Stride { X: 2, Y: 2 }
		Dimensions { K: 512, C: 512, R: 1, S: 1, Y: 28, X: 28 }
	}

	Layer CONV4_1_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 16, C: 512, R: 3, S: 3, Y: 14, X: 14 }
	}

	Layer CONV4_1_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1024, C: 512, R: 1, S: 1, Y: 14, X: 14 }
	}

	Layer CONV4_2_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
	}

	Layer CONV4_2_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 16, C: 512, R: 3, S: 3, Y: 14, X: 14 }
	}

	Layer CONV4_2_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1024, C: 512, R: 1, S: 1, Y: 14, X: 14 }
	}

	Layer CONV4_3_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
	}

	Layer CONV4_3_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 16, C: 512, R: 3, S: 3, Y: 14, X: 14 }
	}

	Layer CONV4_3_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1024, C: 512, R: 1, S: 1, Y: 14, X: 14 }
	}

	Layer CONV4_4_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
	}

	Layer CONV4_4_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 16, C: 512, R: 3, S: 3, Y: 14, X: 14 }
	}

	Layer CONV4_4_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1024, C: 512, R: 1, S: 1, Y: 14, X: 14 }
	}

	Layer CONV4_5_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
	}

	Layer CONV4_5_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 16, C: 512, R: 3, S: 3, Y: 14, X: 14 }
	}

	Layer CONV4_5_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1024, C: 512, R: 1, S: 1, Y: 14, X: 14 }
	}

	Layer CONV4_6_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 512, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
	}

	Layer CONV4_6_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 16, C: 512, R: 3, S: 3, Y: 14, X: 14 }
	}

	Layer CONV4_6_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1024, C: 512, R: 1, S: 1, Y: 14, X: 14 }
	}

	Layer CONV5_1_1 {
		Type: CONV
		Stride { X: 2, Y: 2 }
		Dimensions { K: 1024, C: 1024, R: 1, S: 1, Y: 14, X: 14 }
	}

	Layer CONV5_1_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 32, C: 1024, R: 3, S: 3, Y: 7, X: 7 }
	}

	Layer CONV5_1_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 2048, C: 1024, R: 1, S: 1, Y: 7, X: 7 }
	}

	Layer CONV5_2_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1024, C: 2048, R: 1, S: 1, Y: 7, X: 7 }
	}

	Layer CONV5_2_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 32, C: 1024, R: 3, S: 3, Y: 7, X: 7 }
	}

	Layer CONV5_2_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 2048, C: 1024, R: 1, S: 1, Y: 7, X: 7 }
	}

	Layer CONV5_3_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1024, C: 2048, R: 1, S: 1, Y: 7, X: 7 }
	}

	Layer CONV5_3_2 {
		Type: NGCONV
		Stride { X: 1, Y: 1 }
		Dimensions { G: 32, K: 32, C: 1024, R: 3, S: 3, Y: 7, X: 7 }
	}

	Layer CONV5_3_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 2048, C: 1024, R: 1, S: 1, Y: 7, X: 7 }
	}

	Layer FC1000 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1000, C: 2048, R: 7, S: 7, Y: 7, X: 7 }
	}

}

author: maestro-project

repo: maestro


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

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

%   Copyright 1984-2010 The MathWorks, Inc.
% Last Modified by GUIDE v2.5 30-Sep-2005 11:03:36

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

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

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

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

% Update handles structure
guidata(hObject, handles);

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

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

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

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

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

hMainGui = getappdata(0, 'hMainGui');
fhUpdateAxes = getappdata(hMainGui, 'fhUpdateAxes');

thresh = get(hObject,'Value');
setappdata(hMainGui, 'thresh', thresh);
feval(fhUpdateAxes);

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

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

hMainGui = getappdata(0, 'hMainGui');
thresh   = getappdata(hMainGui, 'thresh');

set(hObject, 'value', thresh)

author: Sable

repo: mcbench-benchmarks


text:
%SETTARGETS Reset targets or soft labels of datafile
%
%     A = SETTARGETS(A,TARGETS)
%
% TARGETS should be a datafile with feature labels
% (FEATLAB) used for the desired LABLIST.

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function int = MITGaussLaguerre( norder, alpha )

%% LAGUERRE_COM computes the abscissa and weights for Gauss-Laguerre
%quadrature.
%
%  Discussion:
%
%    In the simplest case, ALPHA is 0, and we are approximating the
%    integral from 0 to INFINITY of EXP(-X) * F(X).  When this is so,
%    it is easy to modify the rule to approximate the integral from
%    A to INFINITY as well.
%
%    If ALPHA is nonzero, then there is no simple way to extend the
%    rule to approximate the integral from A to INFINITY.  The simplest
%    procedures would be to approximate the integral from 0 to A.
%
%    The integration interval is [ A, +Infinity ) or [ 0, +Infinity ).
%
%    The weight function w(x) = EXP ( - X ) or EXP ( - X ) * X**ALPHA.
%
%    The integral to approximate:
%
%      Integral ( A <= X < +INFINITY ) EXP ( - X ) * F(X) dX
%    or
%      Integral ( 0 <= X < +INFINITY ) EXP ( - X ) * X**ALPHA * F(X) dX
%
%    The quadrature rule:
%
%      EXP ( - A ) * Sum ( 1 <= I <= NORDER ) WEIGHT(I) * F ( A+XTAB(I) )
%    or
%      Sum ( 1 <= I <= NORDER ) WEIGHT(I) * F ( XTAB(I) )
%
%  Modified:
%
%    12 October 2005
%
%  Reference:
%
%    Arthur Stroud and Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966.
%
%  Parameters:
%
%    Input, integer NORDER, the order of the quadrature rule to be
%    computed.
%    NORDER must be at least 1.
%
%    Input, real ALPHA, the exponent of the X factor.
%    Set ALPHA = 0.0 for the simplest rule.
%    ALPHA must be nonnegative.
%
%    Output, real XTAB(NORDER), the Gauss-Laguerre abscissas.
%
%    Output, real WEIGHT(NORDER), the Gauss-Laguerre weights.
%

%
%  Set the recursion coefficients.
%
  for i = 1 : norder
    b(i) = ( alpha + 2 * i - 1 );
  end

  for i = 1 : norder
    c(i) = ( i - 1 ) * ( alpha + i - 1 );
  end

  cc = gamma ( alpha + 1.0 ) * prod ( c(2:norder) );
  int.x = zeros(norder,1);
  int.w = int.x;
  int.n =norder;
  for i = 1 : norder
%
%  Compute an estimate for the root.
%
    if ( i == 1 )

      x = ( 1.0 + alpha ) * ( 3.0+ 0.92 * alpha ) ...
        / ( 1.0 + 2.4 * norder + 1.8 * alpha );

    elseif ( i == 2 )

      x = x + ( 15.0 + 6.25 * alpha ) ...
        / ( 1.0 + 0.9 * alpha + 2.5 * norder );

    else

      r1 = ( 1.0 + 2.55 * ( i - 2 ) ) / ( 1.9 * ( i - 2 ) );

      r2 = 1.26 * ( i - 2 ) * alpha / ( 1.0 + 3.5 * ( i - 2 ) );

      ratio = ( r1 + r2 ) / ( 1.0 + 0.3 * alpha );

      x = x + ratio * ( x - xtab(i-2) );

    end
%
%  Use iteration to find the root.
%
    [ x, dp2, p1 ] = laguerre_root ( x, norder, alpha, b, c );
%
%  Set the abscissa and weight.
%
    int.x(i) = x;
    int.w(i) = ( cc / dp2 ) / p1;

  end
  int.n = norder;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ x, dp2, p1 ] = laguerre_root ( x, norder, alpha, b, c )

%% LAGUERRE_ROOT improves an approximate root of a Laguerre polynomial.
%
%  Modified:
%
%    12 October 2005
%
%  Reference:
%
%    Arthur Stroud and Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966.
%
%  Parameters:
%
%    Input, real X, the approximate root, which
%    should be improved on output.
%
%    Input, integer NORDER, the order of the polynomial to be computed.
%
%    Input, real ALPHA, the exponent of the X factor.
%
%    Input, real B(NORDER), C(NORDER), the recursion coefficients.
%
%    Output, real X, the approximate root, which
%    should be improved on output.
%
%    Output, real DP2, the value of L'(NORDER)(X).
%
%    Output, real P1, the value of L(NORDER-1)(X).
%
  maxstep = 10;

  eps = d_epsilon ( x );

  for i = 1 : maxstep

    [ p2, dp2, p1 ] = laguerre_recur ( x, norder, alpha, b, c );

    d = p2 / dp2;
    x = x - d;

    if ( abs ( d ) <= eps * ( abs ( x ) + 1.0 ) )
      return;
    end

  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ p2, dp2, p1 ] = laguerre_recur ( x, norder, alpha, b, c )

%% LAGUERRE_RECUR finds the value and derivative of a Laguerre
%polynomial.
%
%  Modified:
%
%    12 October 2005
%
%  Reference:
%
%    Arthur Stroud and Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966.
%
%  Parameters:
%
%    Input, real X, the point at which polynomials are evaluated.
%
%    Input, integer NORDER, the order of the polynomial to be computed.
%
%    Input, real ALPHA, the exponent of the X factor in the
%    integrand.
%
%    Input, real B(NORDER), C(NORDER), the recursion
%    coefficients.
%
%    Output, real P2, the value of L(NORDER)(X).
%
%    Output, real DP2, the value of L'(NORDER)(X).
%
%    Output, real P1, the value of L(NORDER-1)(X).
%
  p1 = 1.0;
  dp1 = 0.0;

  p2 = x - alpha - 1.0;
  dp2 = 1.0;

  for i = 2 : norder

    p0 = p1;
    dp0 = dp1;

    p1 = p2;
    dp1 = dp2;

    p2 = ( x - b(i) ) * p1 - c(i) * p0;
    dp2 = ( x - b(i) ) * dp1 + p1 - c(i) * dp0;

  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
function I=getNextCombo(I,n,startVal)
%%GETNEXTCOMBO Return the next combination in lexicographic order given the
%              current combination. If the final combination in the
%              sequence has been reached, then an empty matrix is returned.
%              The first element in the combination is the least
%              significant element for defining the lexicographic order.
%
%INPUTS: I The rX1 or 1Xr current combination of r elements. The next
%          combination in lexicographic order is desired. The first element
%          is the least significant and one begins with I=[0;1;2;...;r-1]
%          is startVal=0 and I=[1;2;3;...;r] if startVal=1.
%        n The number of items from which r items are chosen for
%          combinations. The elements of I can range from 0 to n-1. if
%          startVal=0 or they range from 1 to n if startVal=1.
% startVal This is zero or 1, indicating which value the value at which the
%          elements in I can start. The default if omitted or an empty
%          matrix is passed is 0.
%
%OUTPUTS: I The next combination in the lexicographic sequence, or an empty
%           matrix if the final combination in the lexicographic ordering
%           is provided.
%
%This function can be useful for generating combinations when used in a
%loop. It is more computationally efficient than sequentially unranking the
%combinations. If the final combination is put in, an empty matrix will be
%returned.
%
%The algorithm is from [1], with an option to start at 0 instead of 1.
%
%REFERENCES:
%[1] C. J. Mifsud, "Algorithm 154: Combination in lexicographical order,"
%    Communications of the ACM, vol. 6, no. 3 pp. 103, Mar. 1963.
%    modified to start from zero instead of one.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(startVal))
    startVal=0;
end

if(startVal==0)
    r=length(I);
    if(I(r)<n-1)
        I(r)=I(r)+1;
        return
    else
        for j=r:-1:2
           if(I(j-1)<n-r+j-2)
               I(j-1)=I(j-1)+1;
               for s=j:1:r
                  I(s)=I(j-1)+s-(j-1);
               end
               return;
           end
        end
        I=[];
    end
else%startVal=1
    r=length(I);
    if(I(r)<n)
        I(r)=I(r)+1;
        return
    else
        for j=r:-1:2
           if(I(j-1)<n-r+j-1)
               I(j-1)=I(j-1)+1;
               for s=j:1:r
                  I(s)=I(j-1)+s-(j-1);
               end
               return;
           end
        end
        I=[];
    end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function a = kershaw_llt ( a )

%*****************************************************************************80
%
%% KERSHAW_LLT returns the Cholesky factor of the KERSHAW matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A(4,4), the matrix.
%
  a = [ ...
    1.732050807568877,  0.0,                0.0,                0.0; ...
   -1.154700538379252,  1.290994448735805,  0.0,                0.0; ...
                  0.0, -1.549193338482967,  0.774596669241483,  0.0; ...
    1.154700538379252,  1.032795558988645, -0.516397779494321,  0.577350269189626 ];

  return
end


text:
function mc_install

    if ispc   % A windows version of Matlab
        mex -O -output bbmean winCSource\bbmean.c
        mex -O -output resampres winCsource\resampres.c winCsource\binsgeq.c
        mex -O -output resampsim winCsource\resampsim.c winCsource\binsgeq.c
        mex -O -output resampstr winCsource\resampstr.c winCsource\binsgeq.c
        mex -O -output resampdet winCsource\resampdet.c winCsource\binsgeq.c
    else
        mex -O -output bbmean linuxCsource/bbmean.c
        mex -O -output resampres linuxCsource/resampres.c linuxCsource/binsgeq.c
        mex -O -output resampsim linuxCsource/resampsim.c linuxCsource/binsgeq.c
        mex -O -output resampstr linuxCsource/resampstr.c linuxCsource/binsgeq.c
        mex -O -output resampdet linuxCsource/resampdet.c linuxCsource/binsgeq.c
    end


text:
function x = p05_start ( option, nvar )

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

  return
end


text:
function [Ireg,O_trans,Spacing,M,B,F] = image_registration(Imoving,Istatic,Options)
% This function image_registration is the most easy way to register two
% 2D or 3D images both affine and nonrigidly.
%
% Features:
% - It can be used with images from different type of scans or modalities.
% - It uses both a rigid transform and a nonrigid b-spline grid transform.
% - It uses grid refinement
% - It can be used with images of different sizes.
% - The function will automaticaly detect if the images can be registered
% with the sum of squared pixel distance (SSD), or when mutual information
% must be used as image similarity measure.
%
% Note: Compile the c-files with compile_c_files to allow 3D, and for more
%   more speed in 2D.
%
% [Ireg,Grid,Spacing,M,B,F] = image_registration(Imoving,Istatic,Options);
%
% Inputs,
%   Imoving : The image which will be registerd
%   Istatic : The image on which Imoving will be registered
%   Options : Registration options, see help below
%
% Outputs,
%   Ireg : The registered moving image
%   Grid: The b-spline controlpoints, can be used to transform another
%       image in the same way: I=bspline_transform(Grid,I,Spacing);
%   Spacing: The uniform b-spline knot spacing
%	M : The affine transformation matrix
%   B : The backwards transformation fields  of the pixels in
%       x,y and z direction seen from the  static image to the moving image.
%		in 2D, Bx=B(:,:,1); By=B(:,:,2);
%   F : The (approximated) forward transformation fields of the pixels in
%       x, y and z direction seen from the moving image to the static image.
%		in 2D, Fx=F(:,:,1); Fy=F(:,:,2);
%       (See the function backwards2forwards)
%
% Options,
%   Options.Similarity: Similarity measure (error) used can be set to:
%               sd : Squared pixel distance
%               mi : Normalized (Local) mutual information
%               d, gd, gc, cc, pi, ld : see image_difference.m.
%   Options.Registration:
%               Rigid    : Translation, Rotation
%               Affine   : Translation, Rotation, Shear, Resize
%               NonRigid : B-spline grid based local registration
%               Both     : Nonrigid and Affine (Default)
%   Options.Penalty: Thin sheet of metal smoothness penalty, default in 2D 1e-3 ,
%				default in 3D, 1e-5
%               if set to zero the registration will take a shorter time, but
%               will give a more distorted transformation field.
%   Options.Interpolation: Linear (default) or Cubic, the final result is
%               always cubic interpolated.
%   Options.MaxRef : Maximum number of grid refinements steps (default 2)
%   Options.Grid: Initial B-spline controlpoints grid, if not defined is initalized
%               with an uniform grid. (Used for in example while
%               registering
%               a number of movie frames)
%   Options.Spacing: Spacing of initial B-spline grid in pixels 1x2 [sx sy] or 1x3 [sx sy sz]
%               sx, sy, sz must be powers of 2, to allow grid refinement.
%   Options.MaskMoving: Image which is transformed in the same way as Imoving and
%               is multiplied with the individual pixel errors
%               before calculation of the te total (mean) similarity
%               error. In case of measuring the mutual information similariy
%               the zero mask pixels will be simply discarded.
%   Options.MaskStatic: Also a Mask but is used  for Istatic
%   Options.Verbose: Display Debug information 0,1 or 2
%
% Corresponding points / landmarks Options,
%   Options.Points1: List N x 2 or N x 3 of landmarks x,y(,z) in Imoving image
%   Options.Points2: List N x 2 or N x 3 of landmarks x,y(,z) in Istatic image, in which
%                     every row correspond to the same row with landmarks
%                     in Points1.
%   Options.PStrength: List Nx1 with the error strength used between the
%                     corresponding points, (lower the strenght of the landmarks
%                     if less sure of point correspondence).
%
%
% Example 2D,
%   % Read two greyscale images of Lena
%   Imoving=imread('images/lenag1.png');
%   Istatic=imread('images/lenag2.png');
%
%   % Register the images
%   [Ireg,O_trans,Spacing,M,B,F] = image_registration(Imoving,Istatic);
%
%   % Show the registration result
%   figure,
%   subplot(2,2,1), imshow(Imoving); title('moving image');
%   subplot(2,2,2), imshow(Istatic); title('static image');
%   subplot(2,2,3), imshow(Ireg); title('registerd moving image');
%   % Show also the static image transformed to the moving image
%   Ireg2=movepixels(Istatic,F);
%   subplot(2,2,4), imshow(Ireg2); title('registerd static image');
%
%  % Show the transformation fields
%   figure,
%   subplot(2,2,1), imshow(B(:,:,1),[]); title('Backward Transf. in x direction');
%   subplot(2,2,2), imshow(F(:,:,2),[]); title('Forward Transf. in x direction');
%   subplot(2,2,3), imshow(B(:,:,1),[]); title('Backward Transf. in y direction');
%   subplot(2,2,4), imshow(F(:,:,2),[]); title('Forward Transf. in y direction');
%
% % Calculate strain tensors
%   E = strain(F(:,:,1),F(:,:,2));
% % Show the strain tensors
%   figure,
%   subplot(2,2,1), imshow(E(:,:,1,1),[-1 1]); title('Strain Tensors Exx');
%   subplot(2,2,2), imshow(E(:,:,1,2),[-1 1]); title('Strain Tensors Exy');
%   subplot(2,2,3), imshow(E(:,:,2,1),[-1 1]); title('Strain Tensors Eyx');
%   subplot(2,2,4), imshow(E(:,:,2,2),[-1 1]); title('Strain Tensors Eyy');
%
% Example 2D with Landmarks,
%
% % Read two images with triangles inside
%   Imoving=imread('images/landmarks1.png');
%   Istatic=imread('images/landmarks2.png');
%
% % Load the corresponding (matched) points
%   load('images/landmarks');
%
% % Register the images affine
%   Options = struct('Points1',Points1,'Points2',Points2,'PStrength',PStrength);
%   Ireg = image_registration(Imoving,Istatic,Options);
%
%   % Show the start images
%   figure, imshow(Imoving+Istatic,[]);
%   % Show the registration result
%   figure, imshow(Ireg+Istatic,[]);
%
% Example 2D Prostate,
% % Read two images
%   Imoving=im2double(imread('images/prostate1.png'));
%   Istatic=im2double(imread('images/prostate2.png'));
%
% % Use mutual information
%   Options.Similarity='mi';
% % Set grid smoothness penalty
%   Options.Penalty = 1e-3;
%   Ireg = image_registration(Imoving,Istatic,Options);
%  % Show the registration result
%   figure,
%   subplot(2,2,1), imshow(Imoving); title('moving image');
%   subplot(2,2,2), imshow(Istatic); title('static image');
%   subplot(2,2,3), imshow(Ireg); title('registerd moving image');
%
% Example 3D,
%   % add needed function paths
%   functiondir=which('image_registration.m');
%   addpath([functiondir(1:end-length('image_registration.m')) '/low_level_examples'])
%
%   % Get the volume data
%   [Imoving,Istatic]=get_example_data;
%
%   % Register the images
%   Ireg = image_registration(Imoving,Istatic);
%
%   % Show the results
%   showcs3(Imoving);
%   showcs3(Istatic);
%   showcs3(Ireg);
%
% Function is written by D.Kroon University of Twente (August 2010)

% add all needed function paths
add_function_paths;

% Disable warning
warning('off', 'MATLAB:maxNumCompThreads:Deprecated')

if(size(Imoving,3)>3), IS3D=true; else IS3D=false; end
% Check for presence of needed functions
if(IS3D), check_mexfiles(); end

% Process inputs
defaultoptions=struct('Similarity',[],'Registration','Both','Penalty',1e-3,'MaxRef',2,'Grid',[],'Spacing',[],'MaskMoving',[],'MaskStatic',[],'Verbose',2,'Points1',[],'Points2',[],'PStrength',[],'Interpolation','Linear','Scaling',[1 1]);
if(IS3D), defaultoptions.Penalty=1e-5; defaultoptions.Scaling=[1 1 1]; end
if(~exist('Options','var')), Options=defaultoptions;
else
    tags = fieldnames(defaultoptions);
    for i=1:length(tags),
        if(~isfield(Options,tags{i})), Options.(tags{i})=defaultoptions.(tags{i}); end
    end
    if(length(tags)~=length(fieldnames(Options))),
        warning('image_registration:unknownoption','unknown options found');
    end
end

% Set parameters
type=Options.Similarity;
O_trans=Options.Grid; Spacing=Options.Spacing;
MASKmoving=Options.MaskMoving; MASKstatic=Options.MaskStatic;
Points1=Options.Points1; Points2=Options.Points2; PStrength=Options.PStrength;

% Start time measurement
if(Options.Verbose>0), tic; end

% Convert the input images to double with range 0..1
if(IS3D)
    [Iclass,Imin,Imax,Imoving,Istatic]=images2doublesingle(Imoving,Istatic);
else
    [Iclass,Imin,Imax,Imoving,Istatic]=images2double(Imoving,Istatic);
end

% Resize the moving image to fit the static image
[Istatic,Imoving,MASKmoving]=images2samesize(Istatic,Imoving,MASKmoving,IS3D);

% Detect if the mutual information or pixel distance can be used as
% similarity measure. By comparing the histograms.
if(isempty(type)), type=check_image_modalities(Imoving,Istatic,Options); end

% Register the moving image affine to the static image
if(~strcmpi(Options.Registration(1),'N'))
    M=affine_registration(O_trans,Spacing,Options,Imoving,Istatic,MASKmoving,MASKstatic,type,Points1,Points2,PStrength,IS3D);
else M=[];
end

% Make the initial b-spline registration grid
[O_trans,Spacing,MaxItt]=Make_Initial_Grid(O_trans,Spacing,Options,Imoving,M,IS3D);

% Register the moving image nonrigid to the static image
if(strcmpi(Options.Registration(1),'N')||strcmpi(Options.Registration(1),'B'))
    [O_trans,Spacing]=nonrigid_registration(O_trans,Spacing,Options,Imoving,Istatic,MASKmoving,MASKstatic,type,Points1,Points2,PStrength,MaxItt,IS3D);
end

% Transform the input image with the found optimal grid.
if ( nargout<5 )
    Ireg=bspline_transform(O_trans,Imoving,Spacing,3);
else
    [Ireg,B]=bspline_transform(O_trans,Imoving,Spacing,3);
end

% Make the forward transformation fields from the backwards
if ( nargout>5 ), F=backwards2forwards(B); end

% Convert the double registered image to the class and range of the input images
Ireg=Back2OldRange(Ireg,Iclass,Imin,Imax);

% End time measurement
if(Options.Verbose>0), toc, end

function add_function_paths()
try
    functionname='image_registration.m';
    functiondir=which(functionname);
    functiondir=functiondir(1:end-length(functionname));
    addpath([functiondir '/functions'])
    addpath([functiondir '/functions_affine'])
    addpath([functiondir '/functions_nonrigid'])
catch me
    disp(me.message);
end

function [Istatic,Imoving,MASKmoving]=images2samesize(Istatic,Imoving,MASKmoving,IS3D)
% Resize the moving image to fit the static image
if(sum(size(Istatic)-size(Imoving))~=0)
    % Resize the moving image to fit the static image
    if(IS3D)
        Imoving=imresize3d(Imoving,[],size(Istatic),'cubic');
        if(~isempty(MASKmoving))
            MASKmoving=imresize3d(MASKmoving,[],size(Istatic),'cubic');
        end
    else
        Imoving = imresize(Imoving,[size(Istatic,1) size(Istatic,2)],'bicubic');
        if(~isempty(MASKmoving))
            MASKmoving = imresize(MASKmoving,[size(Istatic,1) size(Istatic,2)],'bicubic');
        end
    end
end

function [Iclass,Imin,Imax,Imoving,Istatic]=images2double(Imoving,Istatic)
% Store the class of the inputs
Iclass=class(Imoving);

% Convert the inputs to double
Imoving=double(Imoving);
Istatic=double(Istatic);
Imin=min(min(Istatic(:)),min(Imoving(:))); Imax=max(max(Istatic(:)),max(Istatic(:)));
Imoving=(Imoving-Imin)/(Imax-Imin);
Istatic=(Istatic-Imin)/(Imax-Imin);

function Ireg=Back2OldRange(Ireg,Iclass,Imin,Imax)
% Back to old image range
Ireg=Ireg*(Imax-Imin)+Imin;

% Set the class of output to input class
if(strcmpi(Iclass,'uint8')), Ireg=uint8(Ireg); end
if(strcmpi(Iclass,'uint16')), Ireg=uint16(Ireg); end
if(strcmpi(Iclass,'uint32')), Ireg=uint32(Ireg); end
if(strcmpi(Iclass,'int8')), Ireg=int8(Ireg); end
if(strcmpi(Iclass,'int16')), Ireg=int16(Ireg); end
if(strcmpi(Iclass,'int32')), Ireg=int32(Ireg); end
if(strcmpi(Iclass,'single')), Ireg=single(Ireg); end

function M=affine_registration(O_trans,Spacing,Options,Imoving,Istatic,MASKmoving,MASKstatic,type,Points1,Points2,PStrength,IS3D)
% Make smooth for fast affine registration
if(IS3D)
    ISmoving=imgaussian(Imoving,2.5,[10 10 10]);
    ISstatic=imgaussian(Istatic,2.5,[10 10 10]);
else
    ISmoving=imfilter(Imoving,fspecial('gaussian',[10 10],2.5));
    ISstatic=imfilter(Istatic,fspecial('gaussian',[10 10],2.5));
end

% Affine register the smoothed images to get the registration parameters
if(strcmpi(Options.Registration(1),'R'))
    if(Options.Verbose>0), disp('Start Rigid registration'); drawnow; end
    if(IS3D)
        x=[0 0 0 0 0 0];
    else
        % Parameter scaling of the Translation and Rotation
        scale=[1 1 1];
        % Set initial affine parameters
        x=[0 0 0];
    end
elseif(strcmpi(Options.Registration(1),'A'))
    if(Options.Verbose>0), disp('Start Affine registration'); drawnow; end
    % Parameter scaling of the Translation, Rotation, Resize and Shear
    if(IS3D)
        x=[0 0 0 0 0 0 1 1 1 0 0 0 0 0 0];
    else
        scale=[1 1 0.01 0.01 0.01 0.01 0.01];
        % Set initial affine parameters
        x=[0 0 0 1 1 0 0];
    end
elseif(strcmpi(Options.Registration(1),'B'))
    if(Options.Verbose>0), disp('Start Affine part of Non-Rigid registration'); drawnow; end
    % Parameter scaling of the Translation, Rotation, Resize and Shear
    if(IS3D)
        x=[0 0 0 0 0 0 1 1 1 0 0 0 0 0 0];
    else
        scale=[1 1 0.01 0.01 0.01 0.01 0.01];
        % Set initial affine parameters
        x=[0 0 0 1 1 0 0];
    end
else
    warning('image_registrations:unknownoption','unknown registration method');
end

if(Options.Interpolation(1)=='L'), interpolation_mode=0; else interpolation_mode=2; end

% Register Affine with 3 scale spaces
for refine_itt=1:3
    if(refine_itt==1)
        if(IS3D)
            ITmoving=imresize3d(ISmoving,0.25);
            ITstatic=imresize3d(ISstatic,0.25);
        else
            ITmoving=imresize(ISmoving,0.25);
            ITstatic=imresize(ISstatic,0.25);
        end
        Points1t=Points1*0.25; Points2t=Points2*0.25; PStrengtht=PStrength;
        if(IS3D)
            if(~isempty(MASKmoving)), ITMASKmoving = imresize3d(MASKmoving,0.25); else ITMASKmoving=[]; end
            if(~isempty(MASKstatic)), ITMASKstatic = imresize3d(MASKstatic,0.25); else ITMASKstatic=[]; end
        else
            if(~isempty(MASKmoving)), ITMASKmoving = imresize(MASKmoving,0.25); else ITMASKmoving=[]; end
            if(~isempty(MASKstatic)), ITMASKstatic = imresize(MASKstatic,0.25); else ITMASKstatic=[]; end
        end
    elseif(refine_itt==2)
        x(1:2)=x(1:2)*2;
        if(IS3D)
            ITmoving=imresize3d(ISmoving,0.5);
            ITstatic=imresize3d(ISstatic,0.5);
        else
            ITmoving=imresize(ISmoving,0.5);
            ITstatic=imresize(ISstatic,0.5);
        end
        Points1t=Points1*0.5; Points2t=Points2*0.5; PStrengtht=PStrength;
        if(IS3D),
            if(~isempty(MASKmoving)), ITMASKmoving = imresize3d(MASKmoving,0.5); else ITMASKmoving=[]; end
            if(~isempty(MASKstatic)), ITMASKstatic = imresize3d(MASKstatic,0.5); else ITMASKstatic=[]; end
        else
            if(~isempty(MASKmoving)), ITMASKmoving = imresize(MASKmoving,0.5); else ITMASKmoving=[]; end
            if(~isempty(MASKstatic)), ITMASKstatic = imresize(MASKstatic,0.5); else ITMASKstatic=[]; end
        end
    elseif(refine_itt==3)
        if(IS3D),
            x(1:3)=x(1:3)*2;
        else
            x(1:2)=x(1:2)*2;
        end
        ITmoving=Imoving;
        ITstatic=Istatic;
        ITMASKmoving = MASKmoving;
        ITMASKstatic = MASKstatic;
        Points1t=Points1; Points2t=Points2; PStrengtht=PStrength;
    end
    % Minimizer parameters
    % Use struct because expanded optimset is part of the Optimization Toolbox.
    if(IS3D)
        optim=struct('GradObj','on','GoalsExactAchieve',1,'Display','off','StoreN',10,'HessUpdate','lbfgs','MaxIter',100,'MaxFunEvals',1000,'TolFun',1e-7,'DiffMinChange',1e-5);
    else
        optim=struct('GradObj','on','GoalsExactAchieve',1,'Display','off','StoreN',10,'HessUpdate','lbfgs','MaxIter',100,'MaxFunEvals',1000,'TolFun',1e-7,'DiffMinChange',1e-3);
    end
    if(Options.Verbose>0), optim.Display='iter'; end
    if(IS3D)
        % Optimal scaling of values, to make the influence of every parameter in the
        % optimizer almost the same:
        [s_Trans, s_Rotation, s_Scale, s_Shear] = affine_parameter_scaling(size(ITmoving));
        if(strcmpi(Options.Registration(1),'R'))
            scale =[s_Trans s_Rotation];
        else
            if(strcmpi(type,'sd')),
                scale =[s_Trans s_Rotation s_Scale s_Shear];
            else
                scale =[s_Trans s_Rotation s_Scale/10 s_Shear/10];
            end
        end
    end

    % Scale the translation, resize and rotation parameters to scaled
    % optimizer values
    x=x./scale;
    % Find the Afffine deformation
    x=fminlbfgs(@(x)affine_registration_error(x,scale,ITmoving,ITstatic,type,O_trans,Spacing,ITMASKmoving,ITMASKstatic,Points1t,Points2t,PStrengtht,interpolation_mode),x,optim);
    % Scale the translation, resize and rotation parameters to the real values
    x=x.*scale;
end

if(strcmpi(Options.Registration(1),'R'))
    % Make the rigid transformation matrix
    if(IS3D)
        M=make_transformation_matrix(x(1:3),x(4:6));
    else
        M=make_transformation_matrix(x(1:2),x(3));
    end
else
    % Make the affine transformation matrix
    if(IS3D)
        M=make_transformation_matrix(x(1:3),x(4:6),x(7:9),x(10:15));
    else
        M=make_transformation_matrix(x(1:2),x(3),x(4:5),x(6:7));
    end
end

function type=check_image_modalities(Imoving,Istatic,Options)
% Detect if the mutual information or pixel distance can be used as
% similarity measure. By comparing the histograms.
Hmoving = hist(Imoving(:),(1/60)*[0.5:60])./numel(Imoving);
Hstatic = hist(Istatic(:),(1/60)*[0.5:60])./numel(Istatic);
if(sum(abs(Hmoving(:)-Hstatic(:)))>0.25),
    type='mi';
    if(Options.Verbose>0), disp('Multi Modalities, Mutual information is used'); drawnow; end
else
    type='sd';
    if(Options.Verbose>0), disp('Same Modalities, Pixel Distance is used'); drawnow; end
end

function [O_trans,Spacing,MaxItt]=Make_Initial_Grid(O_trans,Spacing,Options,Imoving,M,IS3D)
if(isempty(O_trans)),

    if(isempty(Options.Spacing))
        if(IS3D)
            % Calculate max refinements steps
            MaxItt=min(floor(log2(size(Imoving)/4)));

            % set b-spline grid spacing in x,y and z direction
            Spacing=[2^MaxItt 2^MaxItt 2^MaxItt];
        else
            % Calculate max refinements steps
            MaxItt=min(floor(log2([size(Imoving,1) size(Imoving,2)]/4)));

            % set b-spline grid spacing in x and y direction
            Spacing=[2^MaxItt 2^MaxItt];
        end
    else
        % set b-spline grid spacing in x and y direction
        Spacing=round(Options.Spacing);
        t=Spacing; MaxItt=0; while((nnz(mod(t,2))==0)&&(nnz(t<8)==0)), MaxItt=MaxItt+1; t=t/2; end
    end
    % Make the Initial b-spline registration grid
    if(IS3D)
        O_trans=make_init_grid(Spacing,size(Imoving),M);
    else
        if(strcmpi(Options.Registration(1),'N'))
            O_trans=make_init_grid(Spacing,[size(Imoving,1) size(Imoving,2)]);
        else
            O_trans=make_init_grid(Spacing,[size(Imoving,1) size(Imoving,2)],M);
        end
    end
else
    MaxItt=0;
    TestSpacing=Spacing;
    while(mod(TestSpacing,2)==0), TestSpacing=TestSpacing/2; MaxItt=MaxItt+1; end

    if(IS3D)
        % Calculate center of the image
        mean=size(Imoving)/2;
        % Make center of the image coordinates 0,0
        xd=O_trans(:,:,:,1)-mean(1); yd=O_trans(:,:,:,2)-mean(2); zd=O_trans(:,:,:,3)-mean(3);
        % Calculate the rigid transformed coordinates
        O_trans(:,:,:,1) = mean(1) + M(1,1) * xd + M(1,2) *yd + M(1,3) *zd + M(1,4)* 1;
        O_trans(:,:,:,2) = mean(2) + M(2,1) * xd + M(2,2) *yd + M(2,3) *zd + M(2,4)* 1;
        O_trans(:,:,:,3) = mean(3) + M(3,1) * xd + M(3,2) *yd + M(3,3) *zd + M(3,4)* 1;
    else
        if(~strcmpi(Options.Registration(1),'N'))
            % Calculate center of the image
            mean=size(Imoving)/2;
            % Make center of the image coordinates 0,0
            xd=O_trans(:,:,1)-mean(1); yd=O_trans(:,:,2)-mean(2);
            % Calculate the affine transformed coordinates
            O_trans(:,:,1) = mean(1) + M(1,1) * xd + M(1,2) *yd + M(1,3) * 1;
            O_trans(:,:,2) = mean(2) + M(2,1) * xd + M(2,2) *yd + M(2,3) * 1;
        end
    end
end
% Limit refinements steps to user input
if(Options.MaxRef<MaxItt), MaxItt=Options.MaxRef; end

function [O_trans,Spacing]=nonrigid_registration(O_trans,Spacing,Options,Imoving,Istatic,MASKmoving,MASKstatic,type,Points1,Points2,PStrength,MaxItt,IS3D)
% Non-rigid b-spline grid registration
if(Options.Verbose>0), disp('Start non-rigid b-spline grid registration'); drawnow; end

if(IS3D),
    if (Options.Verbose>0), disp(['Current Grid size : ' num2str(size(O_trans,1)) 'x' num2str(size(O_trans,2)) 'x' num2str(size(O_trans,3)) ]); drawnow; end
else
    if (Options.Verbose>0), disp(['Current Grid size : ' num2str(size(O_trans,1)) 'x' num2str(size(O_trans,2))]); drawnow; end
end

% set registration options.
options.type=type;
options.penaltypercentage=Options.Penalty;
options.interpolation=Options.Interpolation;
options.scaling=Options.Scaling;
options.verbose=false;

% Enable forward instead of central gradient incase of error measure is pixel distance
if(strcmpi(type,'sd')), options.centralgrad=false; end

% Make smooth images for fast registration without local minimums
if(IS3D)
    Hsize=round(0.1667*(size(Istatic,1)/size(O_trans,1)+size(Istatic,2)/size(O_trans,2)+size(Istatic,3)/size(O_trans,3)));
    ISmoving=imgaussian(Imoving,Hsize/5,[Hsize Hsize Hsize]);
    ISstatic=imgaussian(Istatic,Hsize/5,[Hsize Hsize Hsize]);
else
    Hsize=round(0.25*(size(Imoving,1)/size(O_trans,1)+size(Istatic,2)/size(O_trans,2)));
    ISmoving=imfilter(Imoving,fspecial('gaussian',[Hsize Hsize],Hsize/5));
    ISstatic=imfilter(Istatic,fspecial('gaussian',[Hsize Hsize],Hsize/5));
end
resize_per=2^(MaxItt-1);

% Modified by Yuan Zhou on 3.7.2017
resize_per = max(resize_per,1);
% Reshape O_trans from a matrix to a vector.
sizes=size(O_trans); O_trans=O_trans(:);

% Resize the mask to the image size used in the registration
if(IS3D)
    if(~isempty(MASKmoving)), MASKmovingsmall=imresize3d(MASKmoving,1/resize_per);  else MASKmovingsmall=[]; end
    if(~isempty(MASKstatic)), MASKstaticsmall=imresize3d(MASKstatic,1/resize_per);  else MASKstaticsmall=[]; end
else
    if(~isempty(MASKmoving)), MASKmovingsmall=imresize(MASKmoving,1/resize_per);  else MASKmovingsmall=[]; end
    if(~isempty(MASKstatic)), MASKstaticsmall=imresize(MASKstatic,1/resize_per);  else MASKstaticsmall=[]; end
end

% Use struct because expanded optimset is part of the Optimization Toolbox.
if(IS3D)
    optim=struct('GradObj','on','GoalsExactAchieve',0,'StoreN',10,'HessUpdate','lbfgs','Display','off','MaxIter',100,'DiffMinChange',0.03,'DiffMaxChange',1,'MaxFunEvals',1000,'TolX',0.05,'TolFun',1e-8);
else
    optim=struct('GradObj','on','GoalsExactAchieve',0,'StoreN',10,'HessUpdate','lbfgs','Display','off','MaxIter',100,'DiffMinChange',0.001,'DiffMaxChange',1,'MaxFunEvals',1000,'TolX',0.005,'TolFun',1e-8);
end
if(Options.Verbose>0), optim.Display='iter'; end

% Start the b-spline nonrigid registration optimizer
if(IS3D)
    ISmoving_small=imresize3d(ISmoving,1/resize_per,[],'linear');
    ISstatic_small=imresize3d(ISstatic,1/resize_per,[],'linear');
else
    ISmoving_small=imresize(ISmoving,1/resize_per);
    ISstatic_small=imresize(ISstatic,1/resize_per);
end
Spacing_small=Spacing/resize_per;
Points1_small=Points1/resize_per;
Points2_small=Points2/resize_per;
PStrength_small=PStrength;
O_trans = resize_per*fminlbfgs(@(x)bspline_registration_gradient(x,sizes,Spacing_small,ISmoving_small,ISstatic_small,options,MASKmovingsmall,MASKstaticsmall,Points1_small,Points2_small,PStrength_small),O_trans/resize_per,optim);

% Reshape O_trans from a vector to a matrix
O_trans=reshape(O_trans,sizes);

for refine_itt=1:MaxItt
    if (Options.Verbose>0), disp('Registration Refinement'); drawnow; end

    % Refine the b-spline grid
    [O_trans,Spacing]=refine_grid(O_trans,Spacing,size(Imoving));

    % Make smooth images for fast registration without local minimums
    if(IS3D)
        Hsize=round(0.1667*(size(ISmoving,1)/size(O_trans,1)+size(ISstatic,2)/size(O_trans,2)+size(ISstatic,3)/size(O_trans,3)));
        ISmoving=imgaussian(Imoving,Hsize/5,[Hsize Hsize Hsize]);
        ISstatic=imgaussian(Istatic,Hsize/5,[Hsize Hsize Hsize]);
    else
        Hsize=round(0.25*(size(ISmoving,1)/size(O_trans,1)+size(ISstatic,2)/size(O_trans,2)));
        ISmoving=imfilter(Imoving,fspecial('gaussian',[Hsize Hsize],Hsize/5));
        ISstatic=imfilter(Istatic,fspecial('gaussian',[Hsize Hsize],Hsize/5));
    end
    resize_per=2^(MaxItt-1-refine_itt);

    % No smoothing in last registration step
    if(IS3D)
        if(refine_itt==MaxItt), ISmoving=Imoving; ISstatic=Istatic; optim.TolX=0.05; resize_per=1; end
    else
        if(refine_itt==MaxItt), ISmoving=Imoving; ISstatic=Istatic; resize_per=1; optim.TolX = 0.03; end
    end

    if(IS3D),
        if (Options.Verbose>0), disp(['Current Grid size : ' num2str(size(O_trans,1)) 'x' num2str(size(O_trans,2)) 'x' num2str(size(O_trans,3)) ]); drawnow; end
    else
        if (Options.Verbose>0), disp(['Current Grid size : ' num2str(size(O_trans,1)) 'x' num2str(size(O_trans,2))]); drawnow; end
    end

    % Reshape O_trans from a matrix to a vector.
    sizes=size(O_trans); O_trans=O_trans(:);

    % Resize the mask to the image size used in the registration
    if(IS3D)
        if(~isempty(MASKmoving)), MASKmovingsmall=imresize3d(MASKmoving,1/resize_per);  else MASKmovingsmall=[]; end
        if(~isempty(MASKstatic)), MASKstaticsmall=imresize3d(MASKstatic,1/resize_per);  else MASKstaticsmall=[]; end
    else
        if(~isempty(MASKmoving)), MASKmovingsmall=imresize(MASKmoving,1/resize_per);  else MASKmovingsmall=[]; end
        if(~isempty(MASKstatic)), MASKstaticsmall=imresize(MASKstatic,1/resize_per);  else MASKstaticsmall=[]; end
    end

    % Start the b-spline nonrigid registration optimizer
    if(IS3D)
        ISmoving_small=imresize3d(ISmoving,1/resize_per,[],'linear');
        ISstatic_small=imresize3d(ISstatic,1/resize_per,[],'linear');
    else
        ISmoving_small=imresize(ISmoving,1/resize_per);
        ISstatic_small=imresize(ISstatic,1/resize_per);
    end
    Spacing_small=Spacing/resize_per;
    Points1_small=Points1/resize_per;
    Points2_small=Points2/resize_per;
    PStrength_small=PStrength;
    O_trans = resize_per*fminlbfgs(@(x)bspline_registration_gradient(x,sizes,Spacing_small,ISmoving_small,ISstatic_small,options,MASKmovingsmall,MASKstaticsmall,Points1_small,Points2_small,PStrength_small),O_trans/resize_per,optim);

    % Reshape O_trans from a vector to a matrix
    O_trans=reshape(O_trans,sizes);
end

function check_mexfiles()
if(exist('bspline_transform_3d_double','file')~=3)
    error('bspline_transform_3d_double mex function not found, compile the c-file');
end
if(exist('bspline_transform_3d_single','file')~=3)
    error('bspline_transform_3d_single mex function not found, compile the c-file');
end

function [Iclass,Imin,Imax,Imoving,Istatic]=images2doublesingle(Imoving,Istatic)
% Store the class of the inputs
Iclass=class(Imoving);

% Convert uint8, uint32 etc. to single.
if(~strcmpi(Iclass,'single')&&~strcmpi(Iclass,'double'))
    Imoving=single(Imoving);
    Istatic=single(Istatic);
end
Imin=min(min(Istatic(:)),min(Imoving(:))); Imax=max(max(Istatic(:)),max(Istatic(:)));
Imoving=(Imoving-Imin)/(Imax-Imin);
Istatic=(Istatic-Imin)/(Imax-Imin);

author: zhouyuanzxcv

repo: Hyperspectral


text:
classdef CEC2017_F22 < PROBLEM
% <single> <real> <constrained>
% CEC'2017 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% G. Wu, R. Mallipeddi, and P. N. Suganthan, Problem definitions and
% evaluation criteria for the CEC 2017 competition on constrained real-
% parameter optimization, National University of Defense Technology, China,
% 2016.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;      % Optimal decision vector
        Mat;	% Rotation matrix
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2017.mat'),'Data');
            obj.O = Data{12}.o;
            obj.M = 1;
            if isempty(obj.D) || obj.D < 30
                obj.D   = 10;
                obj.Mat = Data{12}.M_10;
            elseif obj.D < 50
                obj.D   = 30;
                obj.Mat = Data{12}.M_30;
            elseif obj.D < 100
                obj.D   = 50;
                obj.Mat = Data{12}.M_50;
            else
                obj.D   = 100;
                obj.Mat = Data{12}.M_100;
            end
            obj.lower    = zeros(1,obj.D) - 100;
            obj.upper    = zeros(1,obj.D) + 100;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            Z = Z*obj.Mat';
            PopObj = sum(100*(Z(:,1:end-1).^2-Z(:,2:end)).^2+(Z(:,1:end-1)-1).^2,2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            Z = Z*obj.Mat';
            PopCon(:,1) = sum(Z.^2-10*cos(2*pi*Z)+10,2) - 100;
            PopCon(:,2) = sum(Z,2) - 2*size(Z,2);
            PopCon(:,3) = 5 - sum(Z,2);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function c = reverse_cumsum(v)

c = reverse_vector(cumsum(reverse_vector(v)));

author: kristinbranson

repo: JAABA


text:
function F = spm_Ncdf_jdw(x,u,v)
% Cumulative Distribution Function (CDF) for univariate Normal distributions: J.D.  Williams aproximation
% FORMAT F = spm_Ncdf_jdw(x,u,v)
%
% x - ordinates
% u - mean              [Defaults to 0]
% v - variance  (v>0)   [Defaults to 1]
% F - pdf of N(u,v) at x (Lower tail probability)
%__________________________________________________________________________
%
% spm_Ncdf implements the Cumulative Distribution Function (CDF) for
% the Normal (Gaussian) family of distributions.
%
% References:
%--------------------------------------------------------------------------
% An Approximation to the Probability Integral
% J. D. Williams
% The Annals of Mathematical Statistics, Vol. 17, No. 3. (Sep., 1946), pp.
% 363-365.
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_Ncdf_jdw.m 4836 2012-08-10 15:55:21Z karl $

%-Format arguments
%--------------------------------------------------------------------------
if nargin < 3, v = 1; end
if nargin < 2, u = 0; end

%-Approximate integral
%--------------------------------------------------------------------------
x    = (x - u)./sqrt(abs(v));
F    = sqrt(1 - exp(-(2/pi)*x.^2))/2;
i    = x < 0;
F(i) = -F(i);
F    = F + 1/2;

author: spm


text:
classdef DOC3 < PROBLEM
% <multi> <real> <constrained>
% Benchmark MOP with constraints in both decision and objective spaces

%------------------------------- Reference --------------------------------
% Z. Liu and Y. Wang, Handling constrained multiobjective optimization
% problems with constraints in both the decision and objective spaces. IEEE
% Transactions on Evolutionary Computation, 2019, 23(5): 870-884.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            obj.D = 10;
            obj.lower    = [0 0 0 0 0 0 0 0 0 0.01];
            obj.upper    = [1 1 300 100 200 100 1 100 200 0.03];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values and constraint violations
        function Population = Evaluation(obj,varargin)
            X = varargin{1};
            X = max(min(X,repmat(obj.upper,size(X,1),1)),repmat(obj.lower,size(X,1),1));
            g_temp = -9.* X(:, 6) - 15.* X(:, 9) + 6.* X(:, 2) + 16.* X(:, 3) + 10.* (X(:, 7) + X(:, 8));
            g = (g_temp+400.0551) +1;
            PopObj(:,1) = X(:,1);
            PopObj(:,2) = g.*(1 - (PopObj(:,1))./g);
            % Constraints in objective space
            c(:,1) = max( -(PopObj(:,1).^2 + PopObj(:,2).^2-1), 0);
            c(:,2) = max(-( abs( (-PopObj(:,1) + PopObj(:,2) -0.5)/sqrt(2)) - 0.1/sqrt(2)), 0);
            c(:,3) = max(-( abs( (-PopObj(:,1) + PopObj(:,2) -0)/sqrt(2)) - 0.1/sqrt(2)), 0);
            c(:,4) = max(-( abs( (-PopObj(:,1) + PopObj(:,2) +0.5)/sqrt(2)) - 0.1/sqrt(2)), 0);

            % Constraints in decision space
            c(:,5)  = X(:, 10).* X(:, 4) + 0.02.* X(:, 7) - 0.025.* X(:, 6);
            c(:,6)  = X(:, 10).* X(:, 5) + 0.02.* X(:, 8) - 0.015.* X(:, 9);
            c(:,7)  = abs(X(:, 2) + X(:, 3) - X(:, 4) - X(:, 5)) - 0.0001;
            c(:,8)  = abs(0.03.* X(:, 2) + 0.01.* X(:, 3) - X(:, 10).* (X(:, 4) + X(:, 5))) - 0.0001;
            c(:,9)  = abs(X(:, 4) + X(:, 7) - X(:, 6)) - 0.0001;
            c(:,10) = abs(X(:, 5) + X(:, 8) - X(:, 9)) - 0.0001;
            Population = SOLUTION(X,PopObj,c,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,2);
            R = R./repmat(sqrt(sum(R.^2,2)),1,2);
            R(0.3403<R(:,1) & R(:,1)<0.4782 | 0.6553<R(:,1) & R(:,1)<0.7553 | 0.8782<R(:,1) & R(:,1)<0.9403,:) = [];
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = UniformPoint(100,2);
            R = R./repmat(sqrt(sum(R.^2,2)),1,2);
            R(0.3403<R(:,1) & R(:,1)<0.4782 | 0.6553<R(:,1) & R(:,1)<0.7553 | 0.8782<R(:,1) & R(:,1)<0.9403,:) = nan;
        end
    end
end


text:
function Y = symmetrize(X,grps)
%SYMMETRIZE Symmetrize a tensor X in specified modes.
%
%   Y = symmetrize(X) will symmetrize a tensor X with respect to all
%   modes so that Y is symmetric with respect to any permutation of
%   indices.
%
%   Y = symmetrize(X,MODES) will symmetrize a tensor X with respect to the
%   modes specified by the vector MODES of mode indices. The second
%   argument may alternatively be a cell array of vectors of modes to,
%   e.g., specify that it should be symmetric with respect to mode [1 3] as
%   well as [2 4].
%
%   See also TENSOR, TENSOR/ISSYMMETRIC.
%
%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

%T. Kolda, April 2011
n = ndims(X);
sz = size(X);

% Check that grps exists; if not, create it.
if ~exist('grps','var')
    grps = {1:n};
end

% Check that grps is a cell array.
if ~iscell(grps)
    grps = {grps};
end

if ~isnumeric(grps{1})
    error('MODES must be numeric');
end

% Check tensor dimensions for compatibility with symmetrization
ngrps = length(grps);
for i = 1:ngrps
    dims = grps{i};
    for j = dims(2:end)
        if sz(j) ~= sz(dims(1))
            error('Dimension mismatch for symmetrization');
        end
    end
end

% Check for no overlap in the sets
for i = 1:ngrps
    for j = i+1:ngrps
        if ~isempty(intersect(grps{i},grps{j}))
            error('Cannot haver overlapping symmetries');
        end
    end
end

% Create the combinations for each symmetrized subset
combos = cell(ngrps,1);
for i = 1:ngrps
    combos{i} = perms(grps{i});
end

% Create all the permuations to be averaged
total_perms = prod(cellfun(@length,combos));
sym_perms = repmat(1:n, total_perms, 1);
for i = 1:ngrps
    ntimes = prod(cellfun(@length,combos(1:i-1)));
    ncopies = prod(cellfun(@length,combos(i+1:end)));
    nelems = length(combos{i});

    idx = 1;
    for j = 1:ntimes
        for k = 1:nelems
            for l = 1:ncopies
                sym_perms(idx,grps{i}) = combos{i}(k,:);
                idx = idx + 1;
            end
        end
    end
end

% Create an average tensor
Y = tenzeros(size(X));
for i = 1:total_perms
    Y = Y + permute(X,sym_perms(i,:));
end
Y = Y / total_perms;

% It's not *exactly* symmetric due to oddities in differently ordered
% summations and so on, so let's fix that.
% Idea borrowed from Gergana Bounova:
% http://www.mit.edu/~gerganaa/downloads/matlab/symmetrize.m
for i = 1:total_perms
    Z = permute(Y,sym_perms(i,:));
    Y.data(:) = max(Y.data(:),Z.data(:));
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function ebsd = rotate(ebsd,rot,varargin)
% rotate EBSD
%
% Syntax
%
%   % roate the whoole data set about the z-axis by 90*degree
%   ebsd = rotate(ebsd,10*degree)
%
%   % rotate about the x-axis
%   ebsd = rotate(ebsd,rotation.byAxisAngle(xvector,180*degree))
%
%   % rotate only the spatial data
%   ebsd = rotate(ebsd,180*degree,'keepEuler')
%
%   % rotate about a specific point
%   ebsd = rotate(ebsd,180*degree,'center',[0,0])
%
% Input
%  ebsd - @EBSD
%  angle - double
%  q    - @rotation
%
% Options
%  center - [x,y] center of rotation, default is (0,0)
%
% Flags
%  keepXY    - rotate only the orientation data, i.e. the Euler angles
%  keepEuler - rotate only the spatial data, i.e., the x,y, and z values
%
% Output
%  ebsd - @EBSD

if isa(rot,'double'), rot = rotation.byAxisAngle(zvector,rot); end

% rotate the orientations
if ~check_option(varargin,'keepEuler')
  ebsd.rotations = rotation(rot .* ebsd.rotations);
end

% rotate the spatial data
if ~check_option(varargin,'keepXY')

  center = get_option(varargin,'center',[0,0]);

  ebsd = ebsd - center;

  % store coordinates as vector3d
  V = vector3d(ebsd.prop.x,ebsd.prop.y,0);

  % rotate vertices
  V = rot * V;

  % store back
  ebsd.prop.x = V.x(:);
  ebsd.prop.y = V.y(:);

  ebsd = ebsd + center;

  % rotate the unitcell
  V = vector3d(ebsd.unitCell(:,1),ebsd.unitCell(:,2),0);
  V = rot * V;
  ebsd.unitCell = [V.x(:),V.y(:)];
end

author: mtex-toolbox

repo: mtex

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


text:
classdef (HandleCompatible) OrNode < Node
    % OrNode is a class that represents OR connections in a logical formula
    % For further documentation please have a look at the Node Class.
    % .. Authors
    %     - Thomas Pfau 2016

    properties
    end

    methods
        function res = evaluate(self,assignment, printLevel)
            if ~exist('printLevel','var')
                printLevel = 0;
            end
            res = false;
            for i=1:numel(self.children)
                child = self.children(i);
                if child.evaluate(assignment,printLevel)
                    res = true;
                    break;
                end
            end
            if printLevel >= 1
                fprintf('%s : %i\n',self.toString(),res);
            end
        end

        function cnfNode = convertToCNF(self)
            cnfNode = AndNode();
            for c=1:numel(self.children)
                child = self.children(c);
                CNFChild = child.convertToCNF();
                cnfNode = combineChildren(cnfNode,CNFChild);
                %If the child is again an or node, we need to add all
                %children of that child directly to this node.
            end
        end

        function dnfNode = convertToDNF(self)
            dnfNode = OrNode();
            for c=1:numel(self.children)
                child = self.children(c);
                %If the child is again an or node, we need to add all
                %children of that child directly to this node.
                if isa(child,'OrNode')
                    DNFChild = child.convertToDNF();
                    for cc = 1:numel(DNFChild.children)
                        dnfNode.addChild(DNFChild.children(cc))
                    end
                else
                    dnfNode.addChild(child.convertToDNF());
                end
            end
            %finally, remove all duplicate literal nodes from this node.
            for c = 1:numel(dnfNode.children)
                literals = {};
                childrenToRemove = [];
                childNode = dnfNode.children(c);
                for i = 1 : numel(childNode.children)
                    if isa(childNode.children(i),'LiteralNode')
                        if ~any(~cellfun(@isempty, strfind(literals,childNode.children(i).toString())))
                            literals{end+1} = childNode.children(i).toString();
                        else
                            childrenToRemove(end+1) = i;
                        end
                    end
                end
                childNode.children(childrenToRemove) = [];
            end
        end

        function removeDNFduplicates(self)
            % Assuming this is a DNF head node, removeDNFDuplicates checks
            % all present AND nodes for equality and removes replicates.
            %
            % USAGE:
            %    Node.removeDNFduplicates()
            %
            % OUTPUTS:
            %    Node:    A OrNode with all duplicate And nodes removed.
            %
            i = 1;
            literals = self.getLiterals();
            literals = unique(literals);
            comps = false(numel(self.children),numel(literals));
            for i = 1:numel(self.children)
                comps(i,:) = ismember(literals,self.children(i).getLiterals());
            end
            [~,select] = unique(comps,'rows');
            self.children = self.children(select);
        end

        function res = toString(self,PipeAnd)
            if nargin < 2
                PipeAnd = 0;
            end
            res = '(';
            cstring = '';
            for i=1:numel(self.children)
                child = self.children(i);
                if PipeAnd
                    cstring = [cstring child.toString(PipeAnd) ' | '];
                else
                    cstring = [cstring child.toString(PipeAnd) ' or '];
                end
            end
            if length(cstring) > 1
                if PipeAnd
                    cstring = cstring(1:end-3);
                else
                    cstring = cstring(1:end-4);
                end
            end
            if ~isempty(cstring)
                res = [res cstring ')'];
            else
                res = '';
            end
        end

        function tf = deleteLiteral(self, literalID, keepClauses)
            tf = true;
            if ~exist('keepClauses','var')
                keepClauses = true;
            end
            if ~keepClauses
                % we need to be careful about the nesting
                self.reduce();
            end
            % delete the literals from all non Literal children
            arrayfun(@(x) ~isa(x,'LiteralNode') && x.deleteLiteral(literalID, keepClauses), self.children);
            % now, look for children which are empty, or only contain one
            % element
            toDelete = arrayfun(@(x) (isa(x, 'LiteralNode') && x.contains(literalID) ) || (~isa(x,'LiteralNode') && numel(x.children) <= 1), self.children);
            % and check for one element entries
            mergeChildren = arrayfun(@(x) ~isa(x,'LiteralNode') && numel(x.children) == 1, self.children);
            if any(mergeChildren)
                childsToMerge = self.children(mergeChildren);
                childrenToAdd = OrNode();
                for i = 1:numel(childsToMerge)
                    cchild = childsToMerge(i);
                    childrenToAdd(i) = cchild.children;
                end
                % the following works only on 2017b or newer, but is more
                % efficient.
                % childrenToAdd = arrayfun(@(x) x.children, self.children(mergeChildren));
            end
            self.children(toDelete) = [];
            if exist('childrenToAdd','var')
                for child = 1:numel(childrenToAdd)
                    newChild = childrenToAdd(child);
                    self.children(end+1) = newChild;
                    newChild.parent = self;
                end
            end
            %  fprintf('Removing Literal %s from the following node:\n%s\nLeads to the node:\n%s\n',literalID,originalNodeString,self.toString(1));
        end

        function reduce(self)
            childrenChanged = false;
            mergeNode.children = [];
            for i = 1:numel(self.children)
                cchild = self.children(i);
                cchild.reduce()
                %Check if the child has exactly one child. I
                if numel(cchild.children) == 1
                    %If there is only one child, we can directly add the
                    %child to this node.
                    mergeNode.children = [mergeNode.children,cchild.children];
                    childrenChanged = true;
                elseif isa(cchild,'OrNode')
                    %If its an OR node, we can directly add all children to
                    %this node.
                    mergeNode.children = [mergeNode.children,cchild.children];
                    childrenChanged = true;
                else
                    mergeNode.children = [mergeNode.children,cchild];
                end
            end
            if childrenChanged

                self.children = mergeNode.children;
                for i = 1:numel(self.children)
                    self.children(i).parent = self;
                end
            end
        end

    end

end

author: opencobra

repo: cobratoolbox


text:
function exp = exp_not(exp)
% Negate the given expression.
% Output = Not(Input)
%
% In:
%   Input : An expression
%
% Out:
%   Output : The given expression, negated
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-05-22

if ~exp_beginfun('default','fingerprint_check',0,'fingerprint_create',0,'memoize',0) return; end

exp = ~exp;

exp_endfun;

repo: REST


text:
%% (Internal) Guess ECG signals indexes in a multimodal recording
%
%   [ECG_idx, ECG_header ]= get_ECG_idx_from_header(ECG_header)
%
% Arguments:
%
%      + ECG_header: header info
%
% Output:
%
%      + ECG_idx: indexes of the ECG signals in the multimodal matrix
%
%      + ECG_header: header of a ECG only recording
%
% Example:
%
%
% See also get_PPG_ABP_idx_from_header, ECGwrapper
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function [ECG_idx, ECG_header ]= get_ECG_idx_from_header(ECG_header)

    ECG_idx = [];

    str_12_leads_desc = {'I' 'II' 'III' 'AVR' 'AVL' 'AVF' 'V1' 'V2' 'V3' 'V4' 'V5' 'V6' ... %Standard leads.
                         'MLII' 'MLI' 'MLIII' 'MV1' 'MV2' 'MV3' 'MV4' 'MV5' 'MV6' ... % Modifyied leads.
                         'X' 'Y' 'Z' 'VX' 'VY' 'VZ' ... % Orthogonal leads
                         'E-S' 'A-S' 'A-I' ... % Easy leads
                         };

    % add ECG patterns to be matched in ECG_header.desc
    ECG_pattern_desc = {'ECG' 'BIPOLAR' 'LEAD'};

    if( isfield(ECG_header, 'desc') )
        [~, ~, aux_idx ] = intersect(str_12_leads_desc, upper(strtrim(cellstr(ECG_header.desc))) );
    else
        aux_idx = [];
    end

    if( ECG_header.nsig ~= length(aux_idx) )

        ECG_idx = colvec(aux_idx);

        if( isfield(ECG_header, 'desc') )
            aux_idx = find(any(cell2mat(cellfun(@(b)(cell2mat(cellfun(@(a)(~isempty(strfind(a, b))), rowvec(upper(cellstr(ECG_header.desc))), 'UniformOutput', false))), colvec(ECG_pattern_desc), 'UniformOutput', false)),1));
        else
            aux_idx = [];
        end

        if( ~(isempty(aux_idx) && isempty(ECG_idx)) )
            ECG_idx = [ECG_idx; colvec(aux_idx)];
            Not_ECG_idx = setdiff(1:ECG_header.nsig, ECG_idx);
            if( ~isempty(Not_ECG_idx) )
%                 warning('get_ECG_idx_from_header:SignalsDiscarded', disp_option_enumeration('Some signal/s present are not ECG:', cellstr(ECG_header.desc(Not_ECG_idx,:)) ))
            end
        end

        ECG_idx = sort(ECG_idx);

    else
        % all leads are standard ECG leads
        ECG_idx = 1:ECG_header.nsig;
    end

    if(length(ECG_idx) ~= ECG_header.nsig && nargout > 1)
       %trim the heasig if needed
        ECG_header = trim_ECG_header(ECG_header, ECG_idx);
    end


text:
function perm = ordmmdmex(adjncy) %#ok
% perm = ordmmdmex(adjncy)
%
%   Computes multiple-minimum-degree permutation, for sparse
%   Cholesky. Adjncy is  a sparse symmetric matrix; its diagonal
%   is irrelevant.
%
%   Invokes SPARSPAK-A Release III.
%
% **********  INTERNAL FUNCTION OF CHOLTOOL  **********

%   This file is part of CholTool 1.00
%   Copyright (C) 1998 Jos F. Sturm
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program; if not, write to the Free Software
%   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

%Indicate to the user Matlab cannot find the SeDuMi binaries
sedumi_binary_error();

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [vertices,faces] = mesh_freesurferTRI2matlab(file)

% mesh_freesurferTRI2matlab - Read FreeSurfer tesselation (.tri)
%
% USEAGE: [vertices,faces] = mesh_freesurferTRI2matlab(file)
%
% This function will load an ascii file that contains a one
% line specification of the number of vertices followed
% by rows of vertex points.  It then reads a one line
% specification of the number of faces followed by rows
% of face indices into the vertex rows.  Each vertex row
% contains a vertex index number and 3 x,y,z coordinates.
% Each face row contains a face index and three vertex
% indices.  Vertices in the .tri file are indexed from one
% and those returned are indexed from one.
%
% See also the mesh_freesurfer2matlab function to load
% the tesselations that are created by the mris_convert
% function of freesurfer, which have a different text format
% from those of the BEM .tri files.
%
% The freesurfer tesselations may contain too many faces
% for efficient computations.  If so, try 'reducepatch'.
%
% The returned matrices can be input to the patch command, like so:
%
%    Hpatch = patch('Vertices',vertices,'Faces',faces,...
%                   'EdgeColor',[.8 .8 .8],'FaceColor',[0.9 0.9 0.9]);
%
% This will plot the mesh as a patch object.  See the patch command
% and matlab help for more information on coloring this object.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:57 $

% Licence:  GNU GPL, no implied or express warranties
% History:  03/02 Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fid = fopen(file,'r');

if isequal(fid,-1),
    S=sprintf('Could not open file: "%s"',file);
    error(S);
else

    fprintf('...Reading FreeSurfer Tesselation (.tri)\n');

    tic;

    % Check for comment on first line of file
    frewind(fid); temp = fscanf(fid,'%s',1); frewind(fid);
    if findstr(temp,'#'), temp = fgetl(fid); end

    % Read vertices
    Nvertices = fscanf(fid,'%d',1);
    fprintf('...Reading %d Vertices\n',Nvertices);
    vertices = fscanf(fid,'%f',[4,Nvertices]);
    % remove first row (index) and translate
    vertices = vertices(2:4,:)';

    % Read faces
    Nfaces    = fscanf(fid,'%d',1);
    fprintf('...Reading %d Faces\n',Nfaces);
    faces = fscanf(fid,'%d',[4,Nfaces]);
    % remove first row (index) & translate
    faces = faces(2:4,:)';

    fclose(fid);

    t = toc;
    fprintf('...done (%6.2f sec).\n',t);

end

return

author: PatternRecognition

repo: OpenBMI


text:
%% ------------------------solveAmplitudeFlow.m----------------------------

% Solver for Amplitude Flow as given in Algorithm 1 of the Truncated
% Amplitude Flow (TAF) paper minus the truncations.  The idea was to test the
% amplitude based objective function without removing any outliers. The
% advantage is that we can easily call FASTA. Refer to the userguide for a
% detailed usage of the package.

%  See the script 'testAmplitudeFlow.m' for an example of proper usage of
%  this function.

% PAPER TITLE:
%              Solving Systems of Random Quadratic Equations via Truncated
%              Amplitude Flow.

% ARXIV LINK:
%              https://arxiv.org/pdf/1605.08285.pdf

% INPUTS:
%         A:   Function handle/numerical matrix for data matrix A. The rows
%              of this matrix are the measurement vectors that produce
%              amplitude measurements '\psi'.
%         At:  Function handle/numerical matrix for A transpose.
%         b0:  Observed data vector consisting of amplitude measurements
%              generated from b0 = |A*x|. We assign it to 'psi' to be
%              consistent with the notation in the paper.
%         x0:  The initial vector to be used by any solver.
%        opts: struct consists of the options for the algorithm. For
%              details,see header in solvePhaseRetrieval.m or the User
%              Guide.

% OUPTUT :
%         sol: n x 1 vector. It is the estimated signal.
%        outs: A struct consists of the convergence info. For details,
%              see header in solvePhaseRetrieval.m or the User Guide.

% Note:        When a function handle is used, the value of 'n' (the length
%              of the unknown signal) and 'At' (a function handle for the
%              adjoint of 'A') must be supplied. When 'A' is numeric, the
%              values of 'At' and 'n' are ignored and inferred from the
%              arguments

% DESCRIPTION:
%             The wirtinger flow algorithm uses the squared magnitude
%             objective. This generates a slightly different gradient then
%             the amplitude based objective which is not squared. The TAF
%             paper proposes an amplitude based objective optimization
%             using truncation. This code implements a plain vanilla
%             gradient descent without truncation.

% METHOD:
%         1.) Our implementation uses FASTA, A Forward Backward Splitting
%             Package. FASTA can solve gradient descent schemes provided
%             the objective function expression, its gradient expression,
%             the non smooth term and its proximal operator
%

%         2.) Set the objectve f = @(z) 1/2 * norm(abs(z) - b0)^2.
%             The mask contains the truncated vectors that need to be
%             included.

%         3.) The gradient grad = @(z) (z - b0 .* sign(z)). There is no
%             non smooth term hence g = 0. Conseqently prox = @(z) z
%
%         4.) Send all the above parameters to FASTA which then spits out
%             the estimated signal.

% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

%% -----------------------------START-----------------------------------

function [sol, outs] = solveAmplitudeFlow(A, At, b0, x0, opts)
%     addpath('solvers/linesearch');

    innerOpts = struct;
    innerOpts.maxIters = opts.maxIters;
    innerOpts.maxTime = opts.maxTime;
    innerOpts.tol = opts.tol;
    innerOpts.verbose = opts.verbose;
    innerOpts.recordTimes = opts.recordTimes;
    innerOpts.recordResiduals = opts.recordResiduals;
    innerOpts.recordMeasurementErrors = opts.recordMeasurementErrors;
    innerOpts.recordReconErrors = opts.recordReconErrors;
    innerOpts.xt = opts.xt;

    innerOpts.searchMethod = opts.searchMethod;
    innerOpts.betaChoice = opts.betaChoice;

    [sol, outs] = gradientDescentSolver(A, At, x0, b0, @updateObjective, innerOpts);

    function [f, gradf] = updateObjective(~, ~)
        f = @(z) 1/2 * norm(abs(z) - b0)^2;
        gradf = @(z) (z - b0 .* sign(z));
    end
end

author: tomgoldstein

repo: phasepack-matlab


text:
function f=fun1(p)
f=1.0/(p^2);


text:
function x = discretesample(p, n,varargin)
% Samples from a discrete distribution
%
%   x = discretesample(p, n)
%       independently draws n samples (with replacement) from the
%       distribution specified by p, where p is a probability array
%       whose elements sum to 1.
%
%       Suppose the sample space comprises K distinct objects, then
%       p should be an array with K elements. In the output, x(i) = k
%       means that the k-th object is drawn at the i-th trial.
%
%   Remarks
%   -------
%       - This function is mainly for efficient sampling in non-uniform
%         distribution, which can be either parametric or non-parametric.
%
%       - The function is implemented based on histc, which has been
%         highly optimized by mathworks. The basic idea is to divide
%         the range [0, 1] into K bins, with the length of each bin
%         proportional to the probability mass. And then, n values are
%         drawn from a uniform distribution in [0, 1], and the bins that
%         these values fall into are picked as results.
%
%       - This function can also be employed for continuous distribution
%         in 1D/2D dimensional space, where the distribution can be
%         effectively discretized.
%
%       - This function can also be useful for sampling from distributions
%         which can be considered as weighted sum of "modes".
%         In this type of applications, you can first randomly choose
%         a mode, and then sample from that mode. The process of choosing
%         a mode according to the weights can be accomplished with this
%         function.
%
%   Examples
%   --------
%       % sample from a uniform distribution for K objects.
%       p = ones(1, K) / K;
%       x = discretesample(p, n);
%
%       % sample from a non-uniform distribution given by user
%       x = discretesample([0.6 0.3 0.1], n);
%
%       % sample from a parametric discrete distribution with
%       % probability mass function given by f.
%       p = f(1:K);
%       x = discretesample(p, n);
%

%   Created by Dahua Lin, On Oct 27, 2008
%

%% parse and verify input arguments

assert(isfloat(p), 'discretesample:invalidarg', ...
  'p should be an array with floating-point value type.');

assert(isnumeric(n) && isscalar(n) && n >= 0 && n == fix(n), ...
  'discretesample:invalidarg', ...
  'n should be a nonnegative integer scalar.');

%% main

if numel(p) == 1 && p <= 1 && nargin == 2

  x = ones(1,n);
  return

elseif numel(p) == 1

  if nargin == 2 % without replacement
    if 4*n > p
      rp = randperm(p);
      x = rp(1:n);

    else
      f = zeros(1,p); % flags
      sumf = 0;
      while sumf < n
        f(ceil(p * rand(1,n-sumf))) = 1; % sample w/replacement
        sumf = sum(f); % count how many unique elements so far
      end
      x = find(f > 0);
      x = x(randperm(n));
    end
  else % with replacement
    x = ceil(p * rand(1,n));
  end

  return
end

% process p if necessary
K = numel(p);
if ~isequal(size(p), [1, K])
    p = reshape(p, [1, K]);
end

% construct the bins
edges = [0, cumsum(p)];
s = edges(end);
if abs(s - 1) > eps, edges = edges * (1 / s); end

% draw bins
rv = rand(1, n);
c = histcounts(rv, edges);

% extract samples
xv = find(c);

if numel(xv) == n  % each value is sampled at most once
  x = xv;
else                % some values are sampled more than once
  xc = c(xv);
  d = zeros(1, n);
  dv = [xv(1), diff(xv)];
  dp = [1, 1 + cumsum(xc(1:end-1))];
  d(dp) = dv;
  x = cumsum(d);
end

% randomly permute the sample's order
x = x(randperm(n));