text:
function writeSupplementResultsLatex(ours, geom, ncuts)

gtdir = './figs/supp/gt/';
ourdir = './figs/supp/result/';
ncutdir = './figs/supp/ncut/';
geomdir = './figs/supp/geom/';

nf = numel(ours);

[sval, sind] = sort([ours(:).conservation], 'descend');
ours = ours(sind);
geom = geom(sind);
ncuts = ncuts(sind);

fn = {ours(:).imname};

outname = './results/suppfig.tex';

fid = fopen(outname, 'w');
for f = 1:nf

    if mod(f-1, 6)==0
        fprintf(fid, '\\begin{figure*} \n');
        fprintf(fid, '\\begin{center} \n');
        fprintf(fid, '\\begin{tabular}{cccc}\n');
    end

    fprintf(fid, '\\includegraphics[clip=true,trim=1in 3in 1in 3in, height=1.2in]');
    fprintf(fid, ['{' gtdir]);
    fprintf(fid, [strtok(fn{f}, '.') '_gt.pdf}  &  \\hspace{-.15in} \n ']);
    fprintf(fid, '\\includegraphics[clip=true,trim=1in 3in 1in 3in, height=1.2in]');
    fprintf(fid, ['{' ourdir]);
    fprintf(fid, [strtok(fn{f}, '.') '_res.pdf} & \\hspace{-.15in} \n ']);
    fprintf(fid, '\\includegraphics[clip=true,trim=1in 3in 1in 3in, height=1.2in]');
    fprintf(fid, ['{' geomdir]);
    fprintf(fid, [strtok(fn{f}, '.') '_res.pdf} & \\hspace{-.15in} \n ']);
    fprintf(fid, '\\includegraphics[clip=true,trim=1in 3in 1in 3in, height=1.2in]');
    fprintf(fid, ['{' ncutdir]);
    fprintf(fid, [strtok(fn{f}, '.') '_res.pdf} \\\\ \n']);
    fprintf(fid, ['Gt: (%.2f, %.2f) & Ours: (%.2f,%.2f) & Geom: (%.2f,%.2f) & Ncuts: (%.2f,%.2f) \\\\ \n'], ...
       0, 1, ours(f).efficiency, ours(f).conservation, ...
       geom(f).efficiency, geom(f).conservation, ncuts(f).efficiency, ncuts(f).conservation);

   if (mod(f, 6)==0 || f==nf)
        fprintf(fid, '\\end{tabular} \n');
        fprintf(fid, '\\end{center} \n');
        %fprintf(fid, '\\caption{}\n');
        fprintf(fid, '\\end{figure*}\n');
    end
end
fclose(fid);

author: Cloud-CV

repo: object-proposals


text:
%%Ex 4 Combine looping and branching

sum1 = 0;
sum2 = 0;
N = 9
for k = 1:N
 sum1 = sum1+k;
 if (mod(k,3) == 0)
 sum2 = sum2+k;
 end
end
sum1
sum2

% Output:
%        sum1 = 45
%        sum2 = 18
% Remark: Sum1 = 1+2+3+4+5+6+7+8+9, while sum2 = 3+6+9.

author: TheAlgorithms

repo: MATLAB-Octave


text:
function[varargout]=twodsort(varargin)
%TWODSORT  Distances from data points to nearby grid points.
%
%   [DS,XS,YS]=TWODSORT(X,Y,XO,YO,CUTOFF) returns sorted distances D
%   between data points at locations X,Y and grid points at XO,YO.
%
%   X and Y are arrays of the same size into data point locations. XO and
%   YO are arrays of length M and N, say, specifying the bin center
%   locations of an M x N matrix of grid points, i.e.
%
%      XO= [XO_1 XO_2 ... XO_N]      XO =  [YO_1;
%                                           YO_2;
%                                            ...
%                                           YO_M]
%
%   CUTOFF is the maximum distance to be included in the output arrays.
%
%   The output arrays are M numerical arrays arranged as a length M cell
%   array.  That is, there is one cell per element of Y0. Each numerical
%   array has N columns, i.e., the number of elements of X0, with the
%   number of rows varying between arrays.
%
%   DS gives the distances SQRT((X-XO)^2+(Y-YO)^2) of all data points less
%   than the CUTOFF distance from the (m,n)th grid point, sorted in order
%   of increasing distance.  Entries farther than CUTOFF in all output
%   fields are filled with NaNs.
%
%   XS and YS are corresponding deviations X-XO and Y-YO from the grid
%   point location to each data point.
%
%   The choice to put rows into cell arrays is made because for consistency
%   with SPHERESORT and for convenience in parallelizing POLYSMOOTH.
%   _________________________________________________________________
%
%   Limiting output dimension
%
%   [DS,XS,YS]=TWODSORT(X,Y,XO,YO,[CUTOFF JMAX]), where the fifth input
%   argument is a 2-vector, additionally specifies that number of rows of
%   in each cell of the output will be no larger than JMAX.  This option is
%   useful for the 'fixed population' algorithm in POLYSMOOTH.
%   _________________________________________________________________
%
%   Additional input parameters
%
%   Let's say some additional variables Z1, Z2,...,ZK are given at the data
%   locations X,Y.  Then
%
%   [DS,XS,YS,Z1S,Z2S,...,ZKS]=
%
%                TWODSORT(X,Y,Z1,Z2,...,ZK,XO,YO,CUTOFF);
%
%   also returns the values of these variables.
%
%   Z1S, Z2S,...,ZKS are the same size as the other output arguments, and
%   give the values of Z1, Z2,...,ZK sorted according to distance.
%
%   When there are multiple fields to be mapped, one may instead wish to
%   use the approach described under "One grid, many fields" in POLYSMOOTH.
%   _________________________________________________________________
%
%   See also SPHERESORT, POLYSMOOTH.
%
%   'twodsort --t' runs a test.
%
%   Usage: [ds,xs,ys,indexs]=twodsort(x,y,xo,yo,cutoff);
%          [ds,xs,ys,zs,indexs]=twodsort(x,y,zs,xo,yo,cutoff);
%          [ds,xs,ys,z1s,z2s,...,zNs]=twodsort(x,y,z1,z2,...,zN,xo,yo,cutoff);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2008--2020 J.M. Lilly --- type 'help jlab_license' for details

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

xdata=varargin{1};
ydata=varargin{2};
xo=varargin{end-2};
yo=varargin{end-1};
cutoff=varargin{end};
varargin=varargin(3:end-3);

%In case max # points is not input
Ncutoff=inf;
if length(cutoff)==2
    Ncutoff=cutoff(2);
    cutoff=cutoff(1);
end

xo=xo(:);
yo=yo(:);

K=nargout;
if length(varargin)<K-3
    error('Not enough input arguments.')
end

%K,length(varargin)

if ~aresame(size(xdata),size(ydata))
    error('XDATA and YDATA must be the same size.')
end

for k=1:K
    varargout{k}=cell(length(yo),1);
end
indexo=find(isfinite(xdata)&isfinite(ydata));

if ~isempty(indexo)
    vcolon(xdata,ydata);
    vindex(xdata,ydata,indexo,1);
else
    disp(['No finite data values.']), return
end

%indexall=1:length(xdata);

N=length(xo);
for i=1:length(yo)
    xp=vrep(xdata,N,2)-vrep(xo',length(xdata),1);
    yp=vrep(ydata-yo(i),N,2);
    d=sqrt(xp.^2+yp.^2);
    d(d>cutoff)=nan;
    xp(isnan(d))=nan;
    yp(isnan(d))=nan;

    %size(d)
    if ~allall(~isfinite(d))
        [dsort,sorter]=sort(d,'ascend');
        jj=vrep(1:N,length(xdata),1);
        index=sub2ind(size(d),sorter,jj);

        xp=xp(index);
        yp=yp(index);

        L=min(find(sum(isfinite(dsort),2),1,'last'),Ncutoff);
        varargout{1}{i}=dsort(1:L,:);
        varargout{2}{i}=xp(1:L,:);
        varargout{3}{i}=yp(1:L,:);
        %vsize(d,xp,yp,dsort,sorter)
        for k=4:K
            temp=vrep(varargin{k-3}(indexo),N,2);
            temp=temp(index);
            temp(isnan(dsort))=nan;
            varargout{k}{i}=temp(1:L,:);
        end
    end
end

function[]=twodsort_test

[x,y,z]=peaks;
index=randperm(length(z(:)));
index=index(1:200);

[xdata,ydata,zdata]=vindex(x(:),y(:),z(:),index,1);

%Insert some NANs
xdata(1:7:end)=nan;

xo=(-3:.125:3);
yo=(-3:.125:3);
[xg,yg]=meshgrid(xo,yo);

%[ds,xs,ys,xs2,ys2]=twodsort(xdata,ydata,xdata,ydata,xo,yo,[1 20]);
%vsize(ds,xs,ys,xs2,ys2)

[ds,xs,ys,xs2,ys2]=twodsort(xdata,ydata,xdata,ydata,xo,yo,1);
for i=1:length(xs)
    bool(i)=aresame(xs{i}+vrep(xo,size(xs{i},1),1),xs2{i})&aresame(ys{i}+yo(i),ys2{i});
end
reporttest('TWODSORT consistency',allall(bool));

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function dim_num = p06_dim_num ( )

%*****************************************************************************80
%
%% P06_DIM_NUM returns the spatial dimension for problem p06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer DIM_NUM, the spatial dimension of the
%    dependent variables.
%
  dim_num = 2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = subsasgn( a,s,b )
%SUBSASGN for nested cells
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

if(length(s) > 1)
    posCell = size(s(1).subs,2);
    posCell = s(1).subs{1,posCell};
    posMat = cell2mat(s(2).subs);
    [out,idx] = flattenCellMatrix(a.data);
    out{posCell}(sub2ind(size(out{posCell}),posMat(1),posMat(2),posMat(3))) = b;
    a.data = reconFlatCellMatrix(out,idx);

else
    if(strcmp(s.type,'{}'))
        pos = size(s.subs,2);
        pos = s.subs{1,pos};
        [out,idx] = flattenCellMatrix(a.data);
        out{pos} = b;
        a.data = reconFlatCellMatrix(out,idx);

    else
        pos = size(s.subs,2);
        pos = s.subs{1,pos};
        % a = a.data;
        % meta = a.meta;
        if(isa(b,'TRAFO'))
            if(size(s.subs,2) == 5)
                a.data(pos,:) = b.data; % for kernelImg (due to 5D dataset => cha-cha)
                if iscell(b.meta) ~= 1
                    a.meta = b.meta; %changed for curvelab!
                else
                    a.meta(pos,:) = b.meta;
                end
            else
                a.data(1,pos) = b.data;
                if(~iscell(b.meta) || length(b.meta) > 1)
                    a.meta = b.meta; %changed for curvelab!
                else
                    a.meta(1,pos) = b.meta;
                end
            end
        % a = TRAFO(a,meta);
        else
            [out,idx] = flattenCellMatrix(a.data(1,pos));
            for i=1:length(out)
                out{i}(:) = b;
            end
            a.data(1,pos) = reconFlatCellMatrix(out,idx);
        end
    end
end

end

author: thomaskuestner


text:
function prob_test026 ( )

%*****************************************************************************80
%
%% TEST026 tests BURR_CDF, BURR_CDF_INV, BURR_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST026\n' );
  fprintf ( 1, '  For the Burr PDF:\n' );
  fprintf ( 1, '  BURR_CDF evaluates the CDF;\n' );
  fprintf ( 1, '  BURR_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  BURR_PDF evaluates the PDF;\n' );

  a = 1.0;
  b = 2.0;
  c = 3.0;
  d = 2.0;

  check = burr_check ( a, b, c, d );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST026 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %14f\n', a );
  fprintf ( 1, '  PDF parameter B = %14f\n', b );
  fprintf ( 1, '  PDF parameter C = %14f\n', c );
  fprintf ( 1, '  PDF parameter D = %14f\n', d );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '       X            PDF           CDF            CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ x, seed ] = burr_sample ( a, b, c, d, seed );

    pdf = burr_pdf ( x, a, b, c, d );

    cdf = burr_cdf ( x, a, b, c, d );

    x2 = burr_cdf_inv ( cdf, a, b, c, d );

    fprintf ( 1, ' %14f  %14f  %14f  %14f\n', x, pdf, cdf, x2 );

  end

  return
end


text:
%% visall
% Visualize x-, y-, z-normal slices of a 3D field solution with objects and
% sources.

%%% Syntax
%  visall(scalar3d, [opts])
%  visall(scalar3d, obj_array, [opts])
%  visall(scalar3d, obj_array, src_array, [opts])

%%% Description
% |visall(scalar3d)| vizualizes one of the x-, y-, z-components of the E- or
% H-field stored as an instance of <Scalar3d.html |Scalar3d|>.  The field is
% visualized on x-, y-, z-normal slices, which are shown in 3D as well as in 2D.
%
% |visall(..., obj_array)| visualizes the objects in |obj_array| with the slices
% of a field. The elements of |obj_array| are instances of <EMObject.html
% |EMObject|>.
%
% |vis2d(..., obj_array, src_array)| visualizes the objects and sources in
% |obj_array| and |src_array| with the slices of a field.  The elements of
% |src_array| are instances of <Source.html |Source|>.
%
% The optional argument |opts| is an options structure that controls the
% behavior of visualization.  The allowed fields of |opts| are explained below.
% The default values of the fields are shown in brackets |{}|.
%
% * |opts.withgrid|: |true| or |{false}| to show the grid or not.
% * |opts.withinterp|: |{true}| or |false| to interpolate the field or not.
% * |opts.withpml|: |true| or |{false}| to show the PML regions or not.
% * |opts.withabs|: |true| or |{false}| to show the absolute values of the field
% or not.
% * |opts.cscale|: positive value that controls the color bar range.  Set to
% values less than |1.0| to saturate colors.  The default value is |1.0|.
% * |opts.cmax|: another positive value that controls the color bar range; the
% positive maximum of the color bar range is |cscale * cmax|.  If it is set to
% |NaN|, which is the default, then |cscale * (maximum amplitude of the plotted
% data)| is used for the maximum of the color bar range.
% * |opts.isopaque|: |{true}| or |false| to show opaque slices in 3D or not.
% * |opts.withobjsrc|: |true| or |false| to show the objects and sources or not.
% The default values is |true| when the number of objects is small (<= 20), but
% |false| when the number of objects is large (> 20).  To show only sources
% without objects, provide an empty |obj_array| in |visall()|.
% * |opts.phase|: additional phase angle |phi|.  The field is visualized with an
% additional factor |exp(i*phi)| multiplied.  Not useful if |opts.withabs =
% true|.

%%% Note
% When |opts.isopaque = false| is used, the 3D view shows the slices nicely with
% some transparency.  However, the 2D views show worse-looking images because
% the figure renderer that supports transparency interpolates colors differently
% from the normal renderer.

%%% Example
%   [E, H, obj_array, src_array] = maxwell_run({ARGUMENTS});
%   opts.cscale = 1e-1;
%   opts.wihabs = true;
%   visall(E{Axis.x}, obj_array, src_array, opts);

function visall(scalar3d, varargin)

narginchk(1, 4)
chkarg(istypesizeof(scalar3d, 'Scalar3d'), '"scalar3d" should be instance of Scalar3d.');

iarg = 2;
ivararg = 1;
obj_array = [];
if iarg <= nargin && ~istypesizeof(varargin{ivararg}, 'struct')
	obj_array = varargin{ivararg};
	chkarg(istypesizeof(obj_array, 'EMObject', [1 0]), ...
		'argument %d should be "obj_array" (row vector with EMObject as elements).', iarg);
	iarg = iarg + 1;
	ivararg = ivararg + 1;
end

src_array = [];
if iarg <= nargin && ~istypesizeof(varargin{ivararg}, 'struct')
	src_array = varargin{ivararg};
	chkarg(istypesizeof(src_array, 'Source', [1 0]), ...
		'argument %d should be "src_array" (row vector with Source as elements).', iarg);
	iarg = iarg + 1;
	ivararg = ivararg + 1;
end

no_opts = true;
if iarg <= nargin
	opts = varargin{ivararg};
	chkarg(istypesizeof(opts, 'struct'), 'argument %d should be "opts" (struct).', iarg);
	no_opts = false;
end

if no_opts || ~isfield(opts, 'withgrid')
	opts.withgrid = false;
end
if no_opts || ~isfield(opts, 'withinterp')
	opts.withinterp = true;
end
if no_opts || ~isfield(opts, 'withpml')
	opts.withpml = false;
end
if no_opts || ~isfield(opts, 'withabs')
	opts.withabs = false;
end
if no_opts || ~isfield(opts, 'cscale')
	opts.cscale = 1.0;
end
if no_opts || ~isfield(opts, 'cmax')
	opts.cmax = NaN;
end
if no_opts || ~isfield(opts, 'isopaque')
	opts.isopaque = true;
end
if no_opts || ~isfield(opts, 'withobjsrc')
	opts.withobjsrc = true;
end
if no_opts || ~isfield(opts, 'phase')
	opts.phase = 0;
end

tv = TotalView3d();
tv.scalar3d = scalar3d;
tv.obj_array = obj_array;
tv.src_array = src_array;
tv.withgrid = opts.withgrid;
tv.withinterp = opts.withinterp;
tv.withpml = opts.withpml;
tv.withabs = opts.withabs;
tv.cscale = opts.cscale;
tv.cmax = opts.cmax;
tv.isopaque = opts.isopaque;
tv.withobjsrc = opts.withobjsrc;
tv.phase_angle = opts.phase;

tv.show();

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function erf_test ( )

%*****************************************************************************80
%
%% ERF_TEST tests R4_ERF and R8_ERF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ERF_TEST:\n' );
  fprintf ( 1, '  Test ERF_VALUES, R4_ERF R8_ERF\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X          ERF(X)\n' );
  fprintf ( 1, '                     R4_ERF(X)        Diff\n' );
  fprintf ( 1, '                     R8_ERF(X)        Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx1 ] = erf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_erf ( single ( x ) );
    fx3 = r8_erf ( x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.6g\n', x, fx1 );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end


text:
function TessMat = in_tess_curry( TessFile )
% IN_TESS_CURRY: Read Curry tesselation files

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

% Open tesselation file
fid = fopen(TessFile, 'r');
if fid < 0
    error(['Cannot open file ', TessFile]);
end
% Initialize variables
TessMat = struct('Comment', '', 'Vertices', [], 'Faces', []);
curBlock = [];

% Read file line by line
while 1
    % Read one line
    newLine = fgetl(fid);
    if ~ischar(newLine)
        break;
    end
    % Lines to skip
    if isempty(newLine)
        continue;
    end
    % Identify blocks
    if ~isempty(strfind(newLine, 'LOCATION_LIST START_LIST'))
        curBlock = 'Vertices';
        continue;
    elseif ~isempty(strfind(newLine, 'TRIANGLE_LIST START_LIST'))
        curBlock = 'Faces';
        continue;
    elseif ~isempty(strfind(newLine, 'POINT_DESCRIPTION START_LIST'))
        curBlock = 'Comment';
        continue;
    elseif ~isempty(strfind(newLine, 'END'))
        curBlock = [];
        continue;
    end
    % If no block is open: skip the line
    if isempty(curBlock)
        continue;
    end
    % Numeric values
    if ismember(curBlock, {'Vertices', 'Faces'})
        % Interpret values: 3 values per line
        values = str2num(newLine);
        if (length(values) < 3)
            continue;
        end
        % Add to the list
        TessMat.(curBlock) = [TessMat.(curBlock); values];
    % Text
    else
        TessMat.(curBlock) = [TessMat.(curBlock), newLine];
        % Read only the first line of the comments
        if strcmpi(curBlock, 'Comment')
            curBlock = [];
        end
    end
end
% Close file
fclose(fid);

% Check that some vertices where loaded
if isempty(TessMat.Vertices)
    error('This file does not contain any Curry BEM surface.');
end
% Rebuild tesselation
if isempty(TessMat.Faces)
    % Get the convex envelope of the points
    TessMat.Faces = convhulln(TessMat.Vertices);
else
    % Convert from 0-based to 1-based indices
    TessMat.Faces = TessMat.Faces + 1;
    % Reverse faces order
    TessMat.Faces = TessMat.Faces(:, [1 3 2]);
end

author: brainstorm-tools


text:
function model = rmSearchFit_twoGaussiansPosOnly(model,data,params,wProcess,t)
% rmSearchFit_twoGaussiansPosOnly - wrapper for 'fine' one Gaussian fit
%
% model = rmSearchFit_twoGaussiansPosOnly(model,prediction,data,params);
%
% 2008/01 SOD: split of from rmSearchFit.

% now get original sigmas:
gridSigmas_unique = double(unique(params.analysis.sigmaMajor));
% add upper and lower limit:
expandRange    = double(params.analysis.fmins.expandRange);
gridSigmas = double([0.001.*ones(expandRange,1); ...
    gridSigmas_unique; ...
    params.analysis.maxRF.*ones(expandRange,1)]);

% fminsearch options
searchOptions         = params.analysis.fmins.options;
vethresh              = params.analysis.fmins.vethresh;

% amount of negative fits
nNegFit  = 0;
trends   = t.trends;
t_id     = t.dcid+2;

% convert to double just in case
params.analysis.X = double(params.analysis.X);
params.analysis.Y = double(params.analysis.Y);
params.analysis.allstimimages = double(params.analysis.allstimimages);
data = double(data);

%-----------------------------------
% Go for each voxel
%-----------------------------------
progress = 0;tic;
for ii = 1:numel(wProcess),

    % progress monitor (10 dots)
    if floor(ii./numel(wProcess)*10)>progress,
        % print out estimated time left
        if progress==0,
            esttime = toc.*10;
            if floor(esttime./3600)>0,
                fprintf(1,'[%s]:Estimated processing time: %d voxels: %d hours.\n',...
                    mfilename,numel(wProcess),ceil(esttime./3600));
            else
                fprintf(1,'[%s]:Estimated processing time: %d voxels: %d minutes.\n',...
                    mfilename,numel(wProcess),ceil(esttime./60));
            end;
            fprintf(1,'[%s]:Nonlinear optimization (x,y,sigma):',mfilename);
        end;
        fprintf(1,'.');drawnow;
        progress = progress + 1;
    end;

    % volume index
    vi = wProcess(ii);
    vData = data(:,ii);

    % raw rss value (non-squared) - faster than sum(data(:,vi).^2)
    rawrss     = norm(vData);

    % reset tolFun: Precision of evaluation function. We define RMS
    % improvement relative to the initial raw 'no-fit' data RMS. So, 1
    % means stop if there is less than 1% improvement on the fit:
    searchOptions.TolFun = params.analysis.fmins.options.TolFun.*rawrss;

    % start point from grid fit
    startParams = double([model.x0(vi); model.y0(vi); model.s(vi);...
        model.x02(vi); model.y02(vi); model.s2(vi)]);

    % tight search region [lowerbound upperbound]
    if params.analysis.scaleWithSigmas,
        step = params.analysis.relativeGridStep.*max(startParams([3 6]));
        minstep = params.analysis.maxXY./2./params.analysis.minimumGridSampling;
        step = min(step,minstep);
        maxstep = params.analysis.maxXY./2./params.analysis.maximumGridSampling;
        step = max(step,maxstep);
    else
        step = params.analysis.maxXY./2./params.analysis.maximumGridSampling;
    end;
    boundary.xy    = startParams([1 2 4 5])*[1 1] + ones(4,1)*[-1 1].*(step.*expandRange);

    % interpolated sigmas, so we'll look for the closest one.
    closestvalue = zeros(2,1);
    for n=1:2,
        [tmp,tmp2] = sort(abs(gridSigmas_unique-startParams(3*n)));
        closestvalue(n)       = tmp2(1)+expandRange;
    end
    boundary.sigma     = gridSigmas(closestvalue*[1 1]+[-1 1;-1 1].*expandRange);
    bndParams          = [boundary.xy([1 2],:);boundary.sigma(1,:);...
                          boundary.xy([3 4],:);boundary.sigma(2,:)];

    % actual fitting routine
    if searchOptions.MaxIter>0
        outParams = ...
            fmincon(@(x) rmModelSearchFit_twoGaussiansPosOnly(x,vData,...
            params.analysis.X,...
            params.analysis.Y,...
            params.analysis.allstimimages,trends),...
            startParams,[],[],[],[],bndParams(:,1),bndParams(:,2),...
            @(x) distanceCon(x,startParams,step.*expandRange),searchOptions);
    else
        outParams = startParams;
    end
    %[outParams bndParams(:,1) startParams bndParams(:,2)]

    % make RF, prediction and get rss,b
    Xi = params.analysis.X - outParams(1);   % positive x0 moves center right
    Yi = params.analysis.Y - outParams(2);   % positive y0 moves center up
    rf(:,1) = exp( (Yi.*Yi + Xi.*Xi) ./ (-2.*(outParams(3).^2)) );

    Xi = params.analysis.X - outParams(4);   % positive x0 moves center right
    Yi = params.analysis.Y - outParams(5);   % positive y0 moves center up
    rf(:,2) = exp( (Yi.*Yi + Xi.*Xi) ./ (-2.*(outParams(6).^2)) );

    X = [params.analysis.allstimimages * rf trends];
    b    = pinv(X)*vData;
    rss  = norm(vData-X*b).^2;

    % store results only if the first beta is positive, somehow fmincon
    % outputs negative fits. If the fit is negative keep old (grid) fit.
    if all(b(1:2)>0),
        model.x0(vi)   = outParams(1);
        model.y0(vi)   = outParams(2);
        model.s(vi)    = outParams(3);
        model.x02(vi)  = outParams(4);
        model.y02(vi)  = outParams(5);
        model.s2(vi)   = outParams(6);
        model.rss(vi)  = rss;
        model.b([1 2 t_id],vi)  = b;
    else
        % change the percent variance explained to be just under the
        % current vethresh. So it counts as a 'coarse'-fit but can still be
        % included in later 'fine'-fits
        model.rss(vi)  = (1-max((vethresh-0.01),0)).*model.rawrss(vi);
        nNegFit = nNegFit + 1;
    end
end

% just in case
model.s_major = model.s;
model.s_minor = model.s;

% end time monitor
et  = toc;
if floor(et/3600)>0,
    fprintf(1,'Done [%d hours].\n',ceil(et/3600));
else
    fprintf(1,'Done [%d minutes].\n',ceil(et/60));
end;
fprintf(1,'[%s]:Removed negative fits: %d (%.1f%%).\n',...
    mfilename,nNegFit,nNegFit./numel(wProcess).*100);
return;

%-----------------------------------
% make sure that the pRF can only be moved "step" away from original
% position "startParams" - for the one Gaussian model
function [C, Ceq]=distanceCon(x,startParams,step)
Ceq = [];
dist = x([1 2;4 5])-startParams([1 2;4 5]);
C = sqrt(sum(dist.^2,2)) - step;
return;
%-----------------------------------

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% dim_red_type: dimension reduction type: one of the strings: '', 'pca', 'ica'
% new_dim: relevant only if dim_red_type is not ''
function hglmm(numClusters, is_sampled, dim_red_type, new_dim, numOfEMIter, numOfCpus)

    fv_init;

    if ~exist('numOfEMIter', 'var')
        numOfEMIter = 10;
    end

    if ~exist('numOfCpus', 'var')
        numOfCpus = 4;
    end

    word2vec_sqrt = false;

    mult_factor = LMM_HGLMM_MULT_FACTOR;

    output(1, 'LMM_HGLMM_MULT_FACTOR = %d\n', LMM_HGLMM_MULT_FACTOR);

    vectors_file_name = add_data_dir_base('GoogleNews_vectors_norm.mat');
    output_file_name = add_data_dir_base(sprintf('GoogleNews_norm_hglmm_%d.mat', numClusters));

    if word2vec_sqrt
        vectors_file_name = fname_concat(vectors_file_name, '_sqrt');
        output_file_name = fname_concat(output_file_name, '_sqrt');
    end

    if is_sampled
        vectors_file_name = fname_concat(vectors_file_name, '_sampled');
        output_file_name = fname_concat(output_file_name, '_sampled');
    end

    switch dim_red_type
        case {'pca', 'ica', 'wht'}
            post_fix = sprintf('_%s_%d', dim_red_type, new_dim);
            vectors_file_name = fname_concat(vectors_file_name, post_fix);
            output_file_name = fname_concat(output_file_name, post_fix);
    end

    output(1, 'reading vectors file %s\n', vectors_file_name);
    load(vectors_file_name);
    output(1, 'done\n');

    output(1, 'Calculating HGLMM, %d clusters\n', numClusters);

    % multiply by LMM_HGLMM_MULT_FACTOR=10 due to numerical issue
    % transpose because the rows should be the samples
    vectors_ = mult_factor * vectors_';

    % calculate HGLMM

    output(1, 'numOfCpus = %d\n', numOfCpus);
    output(1, 'numOfEMIter = %d\n', numOfEMIter);

    [m_out, s_out, mu_out, sigma_out, b_out, prior_out, samplesWeightsOut ]=HybridEM(vectors_, numClusters, numOfEMIter, numOfCpus);

    save(output_file_name, 'm_out', 's_out', 'mu_out', 'sigma_out', 'b_out', 'prior_out', 'mult_factor', '-v7.3');

    output(1, 'saved to file: %s\n', output_file_name);

    num_lmm = sum(sum(b_out));
    [rows, cols] = size(b_out);
    lmm_fraction = num_lmm / ( rows*cols );
    output(1, 'num_lmm: %d\n', num_lmm);
    output(1, 'lmm_fraction: %.3f\n', lmm_fraction);

end

author: layumi

repo: Image-Text-Embedding


text:
function [robfitstatus] = canlab_glm_group_levels_run1input(wd, c)
% child process of canlab_glm_group_levels
% (see canlab_glm_README.txt for an overview)
%
% ..
%    Copyright (C) 2013  Luka Ruzic
% ..

load(fullfile(wd,sprintf('env_%04d',c)));
%% PREP

% diaryname = fullfile(wd,sprintf('diary_%04d.log',c));
cd(grpmodeldir)

robfitstatus = 0; %#ok
% grfstatus = 0;

%% robfit
cmd = 'robfit(EXPT, includedcons(c), 0, EXPT.mask)';
% diary(diaryname), fprintf('> %s\n',cmd); diary off
fprintf('> %s\n',cmd);
try
    eval(cmd)
    robfitstatus = 1;
catch exc
    if OPTS.nocatch, cd(STARTINGDIR); rethrow(exc);
    else fprintf('> %s\n',getReport(exc,'extended')); end
    %     else diary(diaryname), fprintf(getReport(exc,'extended')), diary off; end
    robfitstatus = -1;
end

%% grf  % having problems with this running in parallel (display problems?)
% if OPTS.run_grf
%     robdir = filenames(fullfile(grpmodeldir,sprintf('robust%04d',includedcons(c))),'char','absolute');
%     if ~isempty(robdir)
%         announce_string('SIGNIFICANT CLUSTER SIZE ESTIMATION using GRF')
%         try
%             load(fullfile(robdir,'SETUP.mat'));
%
%             cmd = 'sigclext = estimate_cluster_extent(.05, pthresh, SETUP.files);';
%             fprintf('> %s\n',cmd)
%             eval(cmd);
%             close all
%
%             fout = fullfile(robdir,'significant_cluster_extents_grf.txt');
%             dlmwrite(fout,[pthresh' sigclext(:,1)],'precision','%g','delimiter',' '); %#ok
%
%         grfstatus = 1;
%         catch exc
%             if OPTS.nocatch, cd(STARTINGDIR); rethrow(exc)
%             else fprintf('> %s\n',getReport(exc,'extended')); end
%             grfstatus = -1;
%         end
%     end
% end

end

function announce_string(string)

s = sprintf('--  %s  --',string);
l = regexprep(s,'.','-');
fprintf('> \n> \n> \n> %s\n> %s\n> %s\n> \n',l,s,l);

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [grad, grad_W_raw] = B_weighted_average(prev_layers, curr_layer, future_layers)

% raw_weights = prev_layers{1}.a;
input = prev_layers{2}.a;
[dim, nFr] = size(input);

% raw_weights2 = exp(raw_weights);
% weights = raw_weights2 / sum(raw_weights2);

weights = curr_layer.weights;

future_grad = 0;
for i=1:length(future_layers)
    future_grad = future_grad + future_layers{i}.grad;
end

% gradient of the input
% grad_y = sum(future_grad,2);
grad_y = future_grad;
grad = grad_y * weights;

% gradient of the weights
grad_W = input' * grad_y;

if 1
    grad_W_raw = weights .* grad_W' - weights * (weights * grad_W);
else
    grad_W_raw = (diag(weights) - weights' * weights) * grad_W;
    grad_W_raw = grad_W_raw';
end

end

author: singaxiong

repo: SignalGraph


text:
function [rstat]=AddIMUErr(finam, fonam, fenam, nrow, act_row, errdefs, tem_mod, rstat)

%randomize
if (rstat(1)~=0)
    randn('state',rstat);
else
    rstat=randn('state');
end

if isempty(tem_mod) %if no tem_mod is specified, use constant temp=0.
    tem_mod.A=1;
    tem_mod.B=0;
    tem_mod.u=0;
end

%input data
fimu=fopen(finam,'rb');
imu_inp=fread(fimu,[nrow inf],'double');
fclose(fimu);
imu_inp=imu_inp(act_row,:);

ndat=size(imu_inp,2);

%output files
feimu=fopen(fonam,'wb');
ferr=fopen(fenam,'wb');

if (isempty(errdefs))   %do not add error
     out=[0:ndat-1;zeros(1,ndat);imu_inp];
     fwrite(feimu,out,'double');
else    %generate and add errors
    %Generate the system model for the TI imu errors
    [Ati, Bti, Cti, Dti, sPti]=imu_modTI_v000(errdefs);
    nst_ti=size(Ati,1);
    st_ti=sPti*randn(nst_ti,1);

    tem=0;
    tem_dif=0;
    [Atv, Btv, Ctv, sPtv]=imu_modTV_v000(errdefs, tem_dif, tem, 1);
    if (~isempty(Atv))
        nst_tv=size(Atv,1);
        st_tv=sPtv*randn(nst_tv,1);
    else
        nst_tv=0;
    end

    %add errors
    nsen=size(Cti,1);
    if (nst_tv==0) %only ti part
        for in=1:ndat
            imu_out=imu_inp(:,in)+Cti*st_ti+Dti*randn(nsen,1);
            fwrite(ferr,[in-1;st_ti],'double');
            fwrite(feimu,[in-1;tem;imu_out],'double');

            %new error values
            tem=tem_mod.A*tem+tem_mod.B*randn(1)+tem_mod.u;
            st_ti=Ati*st_ti+Bti*randn(nst_ti,1);
        end
    else    %both ti and tv parts are generated
        for in=1:ndat
            imu_out=imu_inp(:,in)+Cti*st_ti+Dti*randn(nsen,1)+Ctv*st_tv;
            fwrite(ferr,[in-1;st_ti;st_tv],'double');
            fwrite(feimu,[in-1;tem;imu_out],'double');

            %new error values
            sr_a=tem_mod.A*tem+tem_mod.B*randn(1)+tem_mod.u;
            tem_dif=sr_a-tem;
            tem=sr_a;
            [Atv, Btv, Ctv]=imu_modTV_v000(errdefs, tem_dif, tem, 0);
            st_tv=Atv*st_tv+Btv*randn(nst_tv,1);
            st_ti=Ati*st_ti+Bti*randn(nst_ti,1);
        end
    end
end
fclose(feimu);
fclose(ferr);

repo: ekfmonoslam


text:
% Codes for CVPR-15 work `Face Alignment by Coarse-to-Fine Shape Searching'
% Any question please contact Shizhan Zhu: zhshzhutah2@gmail.com
% Released on July 25, 2015

function [images,T] = testingsetGeneration(img_root,nameList,bbox,...
    priorModel,priorsInfo,mean_simple_face,target_simple_face)

assert(~xor(nargout >= 3,exist('data','var')));
m = length(nameList);
assert(size(bbox,1)==m);
p = zeros(m,8);
p(:,1:4) = repmat(bbox(:,1),1,4) + ...
    repmat(mean_simple_face(:,1:4),m,1) .* repmat(bbox(:,2)-bbox(:,1),1,4);
p(:,5:8) = repmat(bbox(:,3),1,4) + ...
    repmat(mean_simple_face(:,5:8),m,1) .* repmat(bbox(:,4)-bbox(:,3),1,4);
Tb = getTransToSpecific(p,target_simple_face * priorsInfo.win_size);
ori_images = cell(m,1);
images = cell(m,1);
for i = 1:m
    img = imread([img_root nameList{i}]);
    img = im2uint8(img);
    if size(img,3)>1, img = rgb2gray(img);end;
    ori_images{i} = img;
    images{i} = imtransform(img,Tb(i),'XData',[1 priorsInfo.win_size],...
        'YData',[1,priorsInfo.win_size]);
end;

rotationCenter = priorsInfo.win_size * priorsInfo.rotationCenter;
extractWindow = priorsInfo.win_size * priorsInfo.extractWindow;
assert(extractWindow(4)-extractWindow(3) == extractWindow(2)-extractWindow(1));
extractCellSize = (extractWindow(4)-extractWindow(3)) / priorsInfo.extractCellNum;
predicted = zeros(m,2);
F = zeros(m,priorsInfo.extractCellNum*priorsInfo.extractCellNum*31);
for i = 1:m
    F(i,:) = reshape(...
        vl_hog(single(images{i}(extractWindow(3):extractWindow(4),extractWindow(1):extractWindow(2))),extractCellSize),...
        [1,priorsInfo.extractCellNum*priorsInfo.extractCellNum*31]);
end;
for i = 1:2
    predicted(:,i) = predict(priorModel.B{i},F);
end;

predicted_vote = mean(predicted,2);
ind = find(predicted(:,1) .* predicted(:,2) < 0);
predicted_vote(ind) = 0;
ind = find(predicted(:,1) > priorsInfo.predictedVoteThre & predicted(:,2) > priorsInfo.predictedVoteThre);
predicted_vote(ind) = max(predicted(ind,:),[],2);
ind = find(predicted(:,1) < -priorsInfo.predictedVoteThre & predicted(:,2) < -priorsInfo.predictedVoteThre);
predicted_vote(ind) = min(predicted(ind,:),[],2);

Tr = getTransViaRotateGivenCenter(-predicted_vote,rotationCenter,priorsInfo.rotatorLength);
T = getTransViaMerge(priorsInfo.win_size,Tb,Tr);
for i = 1:m
    images{i} = imtransform(ori_images{i},T(i),'XData',[1 priorsInfo.win_size],...
        'YData',[1 priorsInfo.win_size]);
end;

end

author: zhusz


text:
function r = gen_rand_hash(k,n,m);
persistent buffer
persistent loc
try
    % Preallocate 100 samples everytime and empty that buffer first
    if isempty(buffer) || length(buffer) < loc + n*m
        % Previous approach does not work with parallell toolbox!
        s = rng;
        rng(k);
        buffer = rand(n*m+100,1);
        r = reshape(buffer(1:n*m),n,m);
        loc = 1;
        rng(s);
    end
     r = reshape(buffer(loc:n*m+loc-1),n,m);
     loc = loc + n*m;

catch
    % but rng is not available in all versions...
    s = rand('state');
    rand('state',k)
    r = rand(n,m);
    rand('state',s);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function rgb = labinterp(xi,rgb,xo,varargin)
% Interpolate an sRGB colormap in the Lab colorspace.
%
% Syntax:
%  rgb = labinterp(xi,rgb,xo,<interp1 options>)

M = [...
	+3.2406255,-1.5372080,-0.4986286;...
	-0.9689307,+1.8757561,+0.0415175;...
	+0.0557101,-0.2040211,+1.0569959];
wpt = [0.95047,1,1.08883]; % D65
%
rgb = liRGB2Lab(rgb,M,wpt);
rgb = interp1(xi,rgb,xo,varargin{:});
rgb = liLab2RGB(rgb,M,wpt);
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%labinterp
function rgb = liGammaCor(rgb)
% Gamma correction of sRGB data.
idx = rgb <= 0.0031308;
rgb(idx) = 12.92 * rgb(idx);
rgb(~idx) = real(1.055 * rgb(~idx).^(1/2.4) - 0.055);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%liGammaCor
function rgb = liGammaInv(rgb)
% Inverse gamma correction of sRGB data.
idx = rgb <= 0.04045;
rgb(idx) = rgb(idx) / 12.92;
rgb(~idx) = real(((rgb(~idx) + 0.055) / 1.055).^2.4);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%liGammaInv
function lab = liRGB2Lab(rgb,M,wpt) % Nx3 <- Nx3
% Convert a matrix of sRGB values to Lab.
%
%lab = applycform(rgb,makecform('srgb2lab','AdaptedWhitePoint',wpt));
%
% RGB2XYZ:
xyz = liGammaInv(rgb) / M.';
% Remember to include my license when copying my implementation.
% XYZ2Lab:
xyz = bsxfun(@rdivide,xyz,wpt);
idx = xyz>(6/29)^3;
F = idx.*(xyz.^(1/3)) + ~idx.*(xyz*(29/6)^2/3+4/29);
lab(:,2:3) = bsxfun(@times,[500,200],F(:,1:2)-F(:,2:3));
lab(:,1) = 116*F(:,2) - 16;
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%liRGB2Lab
function rgb = liLab2RGB(lab,M,wpt) % Nx3 <- Nx3
% Convert a matrix of Lab values to sRGB.
%
%rgb = applycform(lab,makecform('lab2srgb','AdaptedWhitePoint',wpt));
%
% Lab2XYZ
tmp = bsxfun(@rdivide,lab(:,[2,1,3]),[500,Inf,-200]);
tmp = bsxfun(@plus,tmp,(lab(:,1)+16)/116);
idx = tmp>(6/29);
tmp = idx.*(tmp.^3) + ~idx.*(3*(6/29)^2*(tmp-4/29));
xyz = bsxfun(@times,tmp,wpt);
% Remember to include my license when copying my implementation.
% XYZ2RGB
rgb = max(0,min(1, liGammaCor(xyz * M.')));
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%cbLab2RGB
% Copyright (c) 2017 Stephen Cobeldick
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%license

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
function [fin_im] =  main_function(or_im)

    % PREPROCESSING
        or_im = pre_pro(or_im);

    % NORMALIZE
        or_im = normal(or_im,1);

    % SEGMENTATION
        blksze = 15; thresh = 0.08;
        [norm_im, mask] = segmentation(or_im, blksze, thresh);

    % DIFFUSION
        norm_im = diffusion(norm_im);

    % ORIENTATION
        [orientim, G_xx, G_yy, G_xy, cos2theta, sin2theta, denom] = orient(norm_im, 1, 5, 6);
    %    showorient(orientim, 20);

    % RELIABILTY
        reliability = reliability_f(G_xx, G_yy, G_xy, cos2theta, sin2theta, denom);

    % FREQUENCY
        [freq] = frequency(norm_im, mask, orientim, 33, 5, 4, 14);

    % GABOR FILTERS
        gabor_im = gabor(freq);

    % MORE FILTERS FOR RIDGE PATTERNS
        new_im = more_filter(norm_im, orientim, gabor_im, 0.5, 0.5);

    % BINARISE
        thres = 0;
        binim = binarise(new_im, thres);

    % Applying reliabilty factor to the final enhanced image
        rel = 0.5;  % reliability factor (above 0.5 is considered to be reliable)
        %binim = binim.*mask.*(reliability>rel);

    fin_im = imcomplement(binim);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [img, fileinfo] = read_hdr(filename)

% load_hdr - loading a radiance RBGE file.
%
%   [img, fileinfo] = read_rle_rgbe(filename);
%
% Written by Lawrence A. Taplin (taplin@cis.rit.edu)
%
% Based loosely on the c-code RGBE implementation written by Bruce Walters
% http://www.graphics.cornell.edu/~bjw/rgbe.html
%
% function to read a run length encoded RGBE picture file and header. This does
% not work if the image is not RLE!!!

do_rescale = 0;

% check if the image is in classical image (eg. TIFF) format
if ~strcmp(filename(end-2:end), 'hdr')
    img = double(imread(filename));
    if do_rescale
        img = rescale(img);
    end
    fileinfo = [];
    return;
end

tic

%open the file
fid = fopen(filename,'r');

if fid<0
    error(['File not found:' filename '.']);
end

%read the magic number
tline = fgetl(fid);
if length(tline)<3 | tline(1:2) ~= '#?'
    error('invalid header');
end
fileinfo.identifier = tline(3:end);

%read the header variables into a structure
tline = fgetl(fid);
while ~isempty(tline)
    %find and set the variable name
    n = strfind(tline,'=');
    if ~isempty(n) % skip stuff I don't understand
        vname = lower(tline(1:n(1)-1));
        vval = tline(n+1:end);
        fileinfo = setfield(fileinfo,vname,tline(n+1:end));
        fprintf('Variable = %s, Value = %s\n',vname, vval);
    end
    %read the next line
    tline = fgetl(fid);
end

%set the resolution variables
tline = fgetl(fid);
fileinfo.Ysign = tline(1);
[fileinfo.height,count,errmsg,nextindex] = sscanf(tline(4:end),'%d',1);
fileinfo.Xsign = tline(nextindex+4);
[fileinfo.width,count,errmsg,nextindex] = sscanf(tline(nextindex+7:end),'%d',1);
fprintf('resolution: %s\n',tline);

%allocate space for the scan line data
img = zeros(fileinfo.height, fileinfo.width, 3);

%read the scanline data
if fileinfo.format == '32-bit_rle_rgbe';
    fprintf('Decoding RLE Data stream\n');
end

%read the remaining data
[data, count] = fread(fid,inf,'uint8');
fclose(fid);

scanline_width = fileinfo.width;
num_scanlines = fileinfo.height;

if ((scanline_width < 8)|(scanline_width > 32767))
    % run length encoding is not allowed so read flat
    img = rgbe2float(reshape(data,fileinfo.width,fileinfo.height,4));
    return;
end

scanline_buffer = repmat(uint8(0),scanline_width,4);
dp = 1; %set the data pointer to the begining of the read data

% read in each successive scanline */
for scanline=1:num_scanlines
    %     if mod(scanline,fix(num_scanlines/100))==0
    %         fprintf('scanline = %d\n',scanline);
    %     end
    %
    if (data(dp) ~= 2) | (data(dp+1) ~= 2)% | (bitget(data(dp+2),8)~=1)
        error('this file is not run length encoded');
    end

    if (bitshift(data(dp+2),8)+data(dp+3))~=scanline_width % -128
        error('wrong scanline width');
    end

    dp = dp+4;

    % read each of the four channels read the scanline into the buffer
    for i=1:4
        ptr = 1;
        while(ptr <= scanline_width)
            if (data(dp) > 128) % a run of the same value
                count = data(dp)-128;
                if ((count == 0)|(count-1 > scanline_width - ptr))
                    warning('bad scanline data');
                end
                scanline_buffer(ptr:ptr+count-1,i) = data(dp+1);
                dp = dp+2;
                ptr = ptr+count;
            else % a non-run
                count = data(dp);
                dp = dp+1;
                if ((count == 0)|(count-1 > scanline_width - ptr))
                    warning('bad scanline data');
                end
                scanline_buffer(ptr:ptr+count-1,i) = data(dp:dp+count-1);
                ptr = ptr+count;
                dp = dp+count;
            end
        end
    end

    % now convert data from buffer into floats
    img(scanline,:,:) = rgbe2float(scanline_buffer);
end

toc

% rescale to 0-1
if do_rescale
    a = min(img(:));
    b = max(img(:));
    img = (img-a)/(b-a);
end

% standard conversion from float pixels to rgbe pixels
% the last dimension is assumed to be color
function [rgbe] = float2rgbe(rgb)
s = size(rgb);
rgb = reshape(rgb,prod(s)/3,3);
rgbe = reshape(repmat(uint8(0),[s(1:end-1),4]),prod(s)/3,4);
v = max(rgb,[],2); %find max rgb
l = find(v>1e-32); %find non zero pixel list
rgbe(l,4) = uint8(round(128.5+log(v)/log(2))); %find E
rgbe(l,1:3) = uint8(rgb(l,1:3)./repmat(2.^(double(rgbe(l,4))-128-8),1,3)); %find rgb multiplier
reshape(rgbe,[s(1:end-1),4]); %reshape back to original dimensions

% standard conversion from rgbe to float pixels */
% note: Ward uses ldexp(col+0.5,exp-(128+8)).  However we wanted pixels */
%       in the range [0,1] to map back into the range [0,1].            */
function [rgb] = rgbe2float(rgbe)
s = size(rgbe);
rgbe = reshape(rgbe,prod(s)/4,4);
rgb = zeros(prod(s)/4,3);
l = find(rgbe(:,4)>0); %nonzero pixel list
rgb(l,:) = double(rgbe(l,1:3)).*repmat(2.^(double(rgbe(l,4))-128-8),1,3);
rgb = reshape(rgb,[s(1:end-1),3]);

author: gpeyre

repo: matlab-toolboxes


text:
function [llcfeat, x, y, wid, hgt] = llc_feature(feature, img, c)
%
% Copyright Aditya Khosla http://mit.edu/khosla
%
% Please cite this paper if you use this code in your publication:
%   A. Khosla, J. Xiao, A. Torralba, A. Oliva
%   Memorability of Image Regions
%   Advances in Neural Information Processing Systems (NIPS) 2012
%

if(~exist('c', 'var'))
	c = conf();
end

[llcfeat, x, y, wid, hgt] = feval(['llc_' feature], img, c);

author: adikhosla

repo: feature-extraction


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Set:
%   Set the specified variable to one.
%   Example: robot = Set(robot, do1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function simulation_open_tool

global configuration robot

robot.tool.tool_open=1;
drawrobot3d(robot, robot.q);

repo: ARTE


text:
function status_out = mlstripcommentsfid(ifid, ofid)
%MLSTRIPCOMMENTSFID Strip comments from code read from file identifier.
%
%   STATUS = MLSTRIPCOMMENTSFID(IFID, OFID) reads MATLAB code, line by line,
%   from the file identifier IFID, strips all MATLAB comments, and writes the
%   result on the output file identifier OFID.  For instance, setting OFID to 1
%   (standard output) will print the result on the screen.  STATUS is 0 if all
%   output was written successfully, and -1 if an error occurred.
%
%   See also MLSTRIPCOMMENTSFILE, MLSTRIPCOMMENTSSTR.

%   Author:      Peter J. Acklam
%   Time-stamp:  2004-03-18 08:23:19 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   % Status is 0 if all output was written successfully, and -1 otherwise.
   status = 0;

   while 1

      % Read input line.
      line = fgetl(ifid);
      if ~ischar(line)
         break;
      end

      % Strip MATLAB comments.
      line = mlstripcommentsstr(line);

      % Write output line.
      count = fprintf(ofid, '%s\n', line);

      % See if an error occurred while writing.
      if count < length(line)
         status = -1;
         break;
      end

   end

   % Only set output variable if it is wanted.
   if nargout
      status_out = status;
   end

author: Sable

repo: mcbench-benchmarks


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%

m = size(X, 1);

for i = 1:m
    deltas = zeros(K, 1);
    rep_examples = repmat(X(i, :), K, 1);
    for j = 1:K
        rep_examples(j, :) = rep_examples(j, :) - centroids(j, :);
        deltas(j, 1) = rep_examples(j, :) * rep_examples(j, :)';
    end
    [min_delta, min_delta_index] = min(deltas);
    idx(i, 1) = min_delta_index;
end

% =============================================================

end

author: rmarquis

repo: coursera-machinelearning


text:
function [a, time_in, time_out, el] = test_resample(a_init, sfreq_out, sfreq_in)
% TEST_RESAMPLE: Test all the methods available in Brainstorm for resamlping.
%
% USAGE:  test_resample(a_init, sfreq_out, sfreq_in)
%         test_resample(a_init)
%         test_resample()
%         [a, time_in, time_out, el] = test_resample();
%
% INPUT:
%     - a_init    : Signal to resample
%                   Default: [cos(t);sin(t)], with t=-pi:.0001:pi
%     - sfreq_out : Output sampling frequency
%                   Default: 4217
%     - sfreq_in  : Initial sampling frequency
%                   Default: 5000
%
% OUTPUT:
%     - a        : Cell-array oo the resampled signals for each method
%     - time_in  : Input time vector (scalar)
%     - time_out : Resampled time vector for each method
%     - el       : Computation time for each method

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

% Define input signal
if (nargin < 1)
    t = -pi:.0001:pi;
    a_init = [cos(t);sin(t)];
end
if (nargin < 3)
    sfreq_in  = 5000;
    sfreq_out = 4217;
end
time_in = (1:size(a_init,2)) / sfreq_in;

% List of available methods
list_methods = {'resample', 'resample-rational', 'resample-cascade', 'interp-decimate-cascade', 'fft-spline'};

% Intialize arrays
a        = cell(1,length(list_methods));
time_out = a;
el       = a;
isFigCreated = 0;
% Loop on all the methods
for i = 1:length(list_methods)
    % Method selection
    method = list_methods{i};
    % Resample
    tic;
    [a{i}, time_out{i}] = process_resample('Compute', a_init, time_in, sfreq_out, method);
    el{i} = toc;
    % Plot initial signal
    if ~isFigCreated
        isFigCreated = 1;
        % Create figure: signal
        hFigSignal = figure('Name', 'Resample: signal', 'NumberTitle', 'off', 'Toolbar', 'figure', 'Units', 'normalized', 'Position', [0 0 1 1]);
        zoom on
        % Plot signal
        hAxesSignal(1) = PlotSignal(hFigSignal, 1, sprintf('Initial signal (%5.3fHz)', sfreq_in), time_in, a_init);
        % Create figure: signal
        hFigSpect(1) = figure('Name', 'Resample: |fft|', 'NumberTitle', 'off', 'Toolbar', 'figure', 'Units', 'normalized', 'Position', [0 0 1 1]);
        zoom on
        % Plot spectrum
        hAxesSpect(1) = PlotSpectrum(hFigSpect, 1, sprintf('Initial signal (%5.3fHz)', sfreq_in), sfreq_in, a_init);
    end
    % Plot resampled signal
    sfreq_out_effective = 1./diff(time_out{i}([1,2]));
    axesTitle = sprintf('%s (%5.3fHz, %3.4fs)', method, sfreq_out_effective, el{i});
    hAxesSignal(i+1) = PlotSignal(hFigSignal, i+1, axesTitle, time_out{i}, a{i});
    % Plot resampled spectrum
    hAxesSpect(i+1) = PlotSpectrum(hFigSpect, i+1, axesTitle, sfreq_out_effective, a{i});
end
% Link all axes
linkaxes(hAxesSignal);
linkaxes(hAxesSpect(2:end));

function hAxes = PlotSignal(hFig, iPlot, axesTitle, t, x)
    hAxes = subplot(2, 4, iPlot, 'Parent', hFig);
    plot(hAxes, t, x);
    title(hAxes, axesTitle);

function hAxes1 = PlotSpectrum(hFig, iPlot, axesTitle, sfreq_in, x)
    % Compute FFT
    ntime = size(x,2);
    nfft = 2^nextpow2(ntime);
    Y = fft(x',nfft)' / ntime;
    f = sfreq_in * linspace(0, 1, nfft);
    % Plot |FFT|
    hAxes1 = subplot(2, 4, iPlot, 'Parent', hFig(1));
    plot(hAxes1, f, 2 * log(max(abs(Y'),1e-5)));
    title(hAxes1, axesTitle);


text:
%filename = 'CantileverSquareSYmmetricMesh';
filename = 'CantileverSquareMedium';
%'CantileverSquareSYmmetricMesh';
%'CantileverSquare';
%filename = 'CantileverSquareSmall';
%'CantileverSquareNew';
%'CantileverSquareNewFine';
%'LshapeTriFine';
%'Lshape';
%'LshapeFine';
%'ArchTriFine';
%'Arch_quad_coarse';
%'Bridge_quad_coarse';
%'BridgeCool_Quadrilateral_Bilinear_Structured_Coarse';
%'Bridge';

ptype = 'MACRO';
initial_case = 'given';
m1 = 0.0101;
m2 = 0.0101;
%cost = {'compliance'};
cost = {'stressNorm'};
%cost = {'stressNorm','compliance'};
%weights = [0.55,0.45];
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'PDE';
constraint_case = 'EQUALITY';

Vfrac_initial = 0.3;
optimality_initial = 1e-5;
constr_initial = 1e-5;

Vfrac_final = 0.3;
optimality_final = 1e-5;
constr_final = 1e-5;

stressNormExponent_initial = 2;
stressNormExponent_final = 32;

optimizer = 'DualNestedInPrimal';
optimizerUnconstrained = 'PROJECTED GRADIENT';

designVariable = 'MicroParams';
ub = 0.989;
lb = 0.011;
homegenizedVariablesComputer = 'ByVademecum';
%
%vademecumFileName = 'SuperEllipseQMax';
%vademecumFileName = 'SuperEllipseQ2';
vademecumFileName = 'SuperEllipseQOptAnalytic';
%
% designVariable = 'Density';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';

kfrac = 2;
nsteps = 20;

plotting = false;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 2000;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% Comparison Of Source Shapes In Different Dimensions
%
% The time-varying pressure signals recorded from a photoacoustic source
% look different depending on the number of dimensions used in the
% simulation. This difference occurs because a point source in 1D
% corresponds to a plane wave in 3D, and a point source in 2D corresponds
% to an infinite line source in 3D. This examples shows the difference
% between the signals recorded in each dimension. It builds on the
% Simulations in One Dimension, Homogeneous Propagation Medium, and
% Simulations in Three Dimensions examples.
%
% author: Bradley Treeby and Ben Cox
% date: 29th January 2011
% last update: 20th October 2011
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% =========================================================================
% SETTINGS
% =========================================================================

% size of the computational grid
Nx = 64;    % number of grid points in the x (row) direction
x = 1e-3;   % size of the domain in the x direction [m]
dx = x/Nx;  % grid point spacing in the x direction [m]

% define the properties of the propagation medium
medium.sound_speed = 1500;      % [m/s]

% size of the initial pressure distribution
source_radius = 2;              % [grid points]

% distance between the centre of the source and the sensor
source_sensor_distance = 10;    % [grid points]

% time array
dt = 2e-9;                      % [s]
t_end = 300e-9;                 % [s]

% computation settings
input_args = {'DataCast', 'single'};

% =========================================================================
% ONE DIMENSIONAL SIMULATION
% =========================================================================

% create the computational grid
kgrid = makeGrid(Nx, dx);

% create the time array
kgrid.t_array = 0:dt:t_end;

% create initial pressure distribution
source.p0 = zeros(Nx, 1);
source.p0(Nx/2 - source_radius:Nx/2 + source_radius) = 1;

% define a single sensor point
sensor.mask = zeros(Nx, 1);
sensor.mask(Nx/2 + source_sensor_distance) = 1;

% run the simulation
sensor_data_1D = kspaceFirstOrder1D(kgrid, medium, source, sensor, input_args{:});

% =========================================================================
% TWO DIMENSIONAL SIMULATION
% =========================================================================

% create the computational grid
kgrid = makeGrid(Nx, dx, Nx, dx);

% create the time array
kgrid.t_array = 0:dt:t_end;

% create initial pressure distribution
source.p0 = makeDisc(Nx, Nx, Nx/2, Nx/2, source_radius);

% define a single sensor point
sensor.mask = zeros(Nx, Nx);
sensor.mask(Nx/2 - source_sensor_distance, Nx/2) = 1;

% run the simulation
sensor_data_2D = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% =========================================================================
% THREE DIMENSIONAL SIMULATION
% =========================================================================

% create the computational grid
kgrid = makeGrid(Nx, dx, Nx, dx, Nx, dx);

% create the time array
kgrid.t_array = 0:dt:t_end;

% create initial pressure distribution
source.p0 = makeBall(Nx, Nx, Nx, Nx/2, Nx/2, Nx/2, source_radius);

% define a single sensor point
sensor.mask = zeros(Nx, Nx, Nx);
sensor.mask(Nx/2 - source_sensor_distance, Nx/2, Nx/2) = 1;

% run the simulation
sensor_data_3D = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});

% =========================================================================
% VISUALISATION
% =========================================================================

figure;
[t_sc, t_scale, t_prefix] = scaleSI(t_end);
plot(kgrid.t_array*t_scale, sensor_data_1D./max(abs(sensor_data_1D)), 'b-');
hold on;
plot(kgrid.t_array*t_scale, sensor_data_2D./max(abs(sensor_data_2D)), 'r-');
plot(kgrid.t_array*t_scale, sensor_data_3D./max(abs(sensor_data_3D)), 'k-');
xlabel(['Time [' t_prefix 's]']);
ylabel('Recorded Pressure [au]');
legend('1D', '2D', '3D');

repo: Aero-matlab


text:
function mv = mv_selectPlotType(type);
%
% mv = mv_selectPlotType(type);
%
% Set the Multi Voxel UI to plot the selected plot
% type.
%
% Default types are:
% 1) Plot the whole time courses for each voxel, after
%    detrending/preprocessing.
% 2) Plot the time courses for each trial and condition,
%    averaged across voxels.
% 3) Plot the time courses for each voxel and condition,
%    averaged across trials.
% 4) Plot the mean amplitudes for each voxel and condition,
%    averaged across trials.
% 5) Plot the mean amplitudes for each trial, voxel, and condition.
% 6) d' histograms [n. y. i.]
% 7) Visualize GLMs for each voxel (w/ slider to page between voxels)
%
% ras, 04/05.
mv = get(gcf,'UserData');
mv.ui.plotType = type;
set(gcf, 'UserData', mv);

% clear previous objects in the figure
old = findobj('Type', 'axes', 'Parent', gcf);
old = [old; findobj('Type', 'uicontrol', 'Parent', gcf)];
old = [old; findobj('Tag', 'GLM Display Panel', 'Parent', gcf)];
delete(old);

% hide GLM control panel if we're not visualizing this
if type ~= 7
    mrvPanelToggle(mv.ui.glmPanel, 'off');
end

% make only the selected option checked
if checkfields(mv, 'ui', 'plotHandles')
    set(mv.ui.plotHandles, 'Checked', 'off')
    set(mv.ui.plotHandles(type), 'Checked', 'on');
end

% multiVoxelUI;

return


text:
function varargout = ellipseAsPolygon(ellipse, N)
%ELLIPSEASPOLYGON Convert an ellipse into a series of points.
%
%   Deprecated, use ellipseToPolygon instead.
%
%   P = ellipseAsPolygon(ELL, N);
%   converts ELL given as [x0 y0 a b] or [x0 y0 a b theta] into a polygon
%   with N edges. The result P is (N+1)-by-2 array containing coordinates
%   of the N+1 vertices of the polygon.
%   The resulting polygon is closed, i.e. the last point is the same as the
%   first one.
%
%   P = ellipseAsPolygon(ELL);
%   Use a default number of edges equal to 72. This result in one piont for
%   each 5 degrees.
%
%   [X Y] = ellipseAsPolygon(...);
%   Return the coordinates o fvertices in two separate arrays.
%
%   See also:
%   ellipses2d, circleAsPolygon, rectAsPolygon, drawEllipse

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2005-04-06
% Copyright 2005 INRA - TPV URPOI - BIA IMASTE

warning('matGeom:deprecated', ...
    'function "ellipseAsCurve" is deprecated, use "ellipseToPolygon" instead');

% format output
if nargout <= 1
    varargout = {ellipseToPolygon(ellipse, N)};
else
    [x, y] = ellipseToPolygon(ellipse, N);
    varargout = {x, y};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% PartitionTable.m by David Terr, Raytheon, 6-7-04

% Given nonnegative integer n, compute the unrestricted partition function p(k) for k<=n.
function pt = partitiontable(n)

if n==0
    pt = 1;
    return;
end

pt = zeros( n+1, 2 );
pt( n+1, 1 ) = n;

% Precompute tables of (-1)^k and k*(3k +/- 1)/2.
s = sqrt( 24*n + 1 );
kp = floor( ( s - 1 ) / 6 );
km = floor( ( s + 1 ) / 6 );
tp = zeros( kp );
tm = zeros( km );

for k = 1:kp
    tp( k ) = k*(3*k+1)/2;
end

for k = 1:km
    tm( k ) = k*(3*k-1)/2;
end

% Compute p(m) for all m <= n.
pt( 1, 2 ) = 1;

for m=1:n
    pt(m,1) = m-1;
    s = sqrt( 24*m + 1 );
    kp = floor( ( s - 1 ) / 6 );
    km = floor( ( s + 1 ) / 6 );

    for k = 1:kp
        pt( m+1, 2 ) = pt( m+1, 2 ) + (-1)^(k+1) * pt( m + 1 - tp(k), 2 );
    end

    for k = 1:km
        pt( m+1, 2 ) = pt( m+1, 2 ) + (-1)^(k+1) * pt( m + 1 - tm(k), 2 );
    end
end


text:
%the object be controlled
function [sys,x0,str,ts]=ctrl5_2_3obj(t,x,u,flag)
switch flag,
    case 0,
        [sys,x0,str,ts]=mdlInit();
    case 1,
        sys=mdlDer(t,x,u);
    case 3,
        sys=mdlOutput(t,x,u);
    otherwise,
        sys=[];
end;
function [sys,x0,str,ts]= mdlInit()
size=simsizes;
size.NumContStates=2;
size.NumDiscStates=0;
size.NumOutputs=2;
size.NumInputs=1;
size.DirFeedthrough=0;
size.NumSampleTimes=1;
sys=simsizes(size);
x0=[0];
str=[];
ts=[0,0];
function sys=mdlDer(t,x,u)
global b0;
a=1;
sys(1)=x(2);
sys(2)=a*sign(sin(10*t))+u*b0;
function sys=mdlOutput(t,x,u)
sys(1)=x(1);
sys(2)=sign(sin(2*t));

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
% 2021 EDIT: This function is now deprecated. The current version of the
% shaped PRESS simulation is sim_press_shaped.m, which employs
% coherence selection rather than phase cycling to null undesired signal.

% sim_press_shaped_phCyc.m
% Robin Simpson and Jamie Near, 2014.
%
% USAGE:
% out = sim_press_shaped_phCyc(n,sw,Bfield,linewidth,sys,tau1,tau2,RF,tp,dx,dy,Gx,Gy,phCyc1,phCyc2,flipAngle,centreFreq)
%
% DESCRIPTION:
% This function simulates the PRESS experiment.  The excitation is
% simulated as an instantaneous rotation, and the refocusing pulse is
% simulated as a shaped rotation.
%
% This code enables the choice of the phase of the refocusing pulses.  This
% enables phase cycling of the refocusing pulses by repeating simulations
% with different editing pulse phases, which is necessary to remove phase
% artefacts from the editing pulses.  A four step phase cycling scheme is typically
% sufficient, where both refocusing pulses are phase cycled by 0 and 90 degrees, and
% the phase are combined in the following way:
%
% signal = ([0 90] - [0 0]) + ([90 0] - [90 90]);
%
% where, in [X Y], X is the phase of the first refocusing pulse and Y is
% the phase of the second refocusing pulse
%
% Finally, this code simulates the spectrum at a given point in space (x,y),
% given the values of the slice selection gradients (Gx, and Gy).  The pulse
% waveform is assumed to be the same for both refocusing pulses.  In order
% to fully simulate the MEGA-PRESS experiment, you have to run this
% simulation many times at various points in space (x,y), and then add
% together the resulting spectra.
%
% INPUTS:
% n         = number of points in fid/spectrum
% sw        = desired spectral width in [Hz]
% Bfield    = main magnetic field strength in [T]
% linewidth = linewidth in [Hz]
% sys       = spin system definition structure
% tau1      = echo time 1 in [ms].
% tau2      = echo time 2 in [ms].
% RF        = RF pulse definition structure for refoc pulses (obtain using 'io_loadRFwaveform.m')
% tp        = RF pulse duration in [ms]
% dx        = position offset in x-direction (corresponding to first refocusing pulse) [cm]
% dy        = position offset in y-direction (corresponding to second refocusing pulse) [cm]
% Gx        = gradient strength for first selective refocusing pulse [G/cm]
% Gy        = gradient strength for second selective refocusing pulse [G/cm]
% phCycl    = initial phase of the first refocusing pulse in [degrees];
% phCycl2   = initial phase of the second refocusing pulse in [degrees];
% flipAngle = flip angle of refocusing pulses [degrees] (Optional.  Default = 180 deg)
% centreFreq= centre frequency of the spectrum in [ppm] (Optional.  Default = 2.3)
%
% OUTPUTS:
% out       = simulated spectrum, in FID-A structure format, using PRESS
%             sequence.

function out = sim_press_shaped_phCyc(n,sw,Bfield,linewidth,sys,tau1,tau2,RF,tp,dx,dy,Gx,Gy,phCyc1,phCyc2,flipAngle,centreFreq)

if nargin<17
    centreFreq=2.3;
    if nargin<16
        flipAngle=180;
    end
end

if tau1<tp/1000
    error('ERROR:  Echo-time 1 cannot be less than duration of refocusing pulse! ABORTING!!');
end
if tau2<tp/1000
    error('ERROR:  Echo-time 2 cannot be less than duration of refocusing pulse! ABORTING!!');
end

%Set water to centre
%centreFreq=2.3;
for k=1:length(sys)
    sys(k).shifts=sys(k).shifts-centreFreq;
end

%Calculate Hamiltonian matrices and starting density matrix.
[H,d]=sim_Hamiltonian(sys,Bfield);

%Calculate new delays by subtracting the pulse duration from tau1 and tau2;
delays=zeros(2);
delays(1)=tau1-tp;
delays(2)=tau2-tp;
if sum(delays<0)
    error(['ERROR! The following taus are too short: ' num2str(find(delays<0)) '.']);
end

%BEGIN PULSE SEQUENCE************
d=sim_excite(d,H,'x');                                    %EXCITE
d=sim_evolve(d,H,delays(1)/2000);                            %Evolve by delays(1)/2
d=sim_shapedRF(d,H,RF,tp,flipAngle,90+phCyc1,dx,Gx);          %1st shaped 180 degree refocusing pulse
d=sim_evolve(d,H,(delays(1)+delays(2))/2000);                     %Evolve by (delays(1)+delays(2))/2
d=sim_shapedRF(d,H,RF,tp,flipAngle,90+phCyc2,dy,Gy);          %2nd shaped 180 degree refocusing pulse
d=sim_evolve(d,H,delays(2)/2000);                            %Evolve by delays(2)/2
[out,dout]=sim_readout(d,H,n,sw,linewidth,90);      %Readout along y (90 degree phase);
%END PULSE SEQUENCE**************

%Correct the ppm scale:
out.ppm=out.ppm-(4.65-centreFreq);

%Fill in structure header fields:
out.seq='press';
out.te=tau1+tau2;
out.sim='shaped';

%Additional fields for compatibility with FID-A processing tools.
out.sz=size(out.specs);
out.date=date;
out.dims.t=1;
out.dims.coils=0;
out.dims.averages=0;
out.dims.subSpecs=0;
out.dims.extras=0;
out.averages=1;
out.rawAverages=1;
out.subspecs=1;
out.rawSubspecs=1;
out.flags.writtentostruct=1;
out.flags.gotparams=1;
out.flags.leftshifted=0;
out.flags.filtered=0;
out.flags.zeropadded=0;
out.flags.freqcorrected=0;
out.flags.phasecorrected=0;
out.flags.averaged=1;
out.flags.addedrcvrs=1;
out.flags.subtracted=1;
out.flags.writtentotext=0;
out.flags.downsampled=0;
out.flags.isFourSteps=0;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function sqrtPd = getSqrtPrecon(problem, x, d, storedb, key)
% Applies the square root of the Hessian preconditioner at x along d.
%
% function sqrtPd = getSqrtPrecon(problem, x, d)
% function sqrtPd = getSqrtPrecon(problem, x, d, storedb)
% function sqrtPd = getSqrtPrecon(problem, x, d, storedb, key)
%
% Returns as sqrtPd the result of applying the square root of the Hessian
% preconditioner to the tangent vector d at point x. The preconditioner is
% supposed to be a symmetric, positive definite approximation of the
% inverse of the Hessian. Its square root must thus be symmetric and
% positive definite itself.
%
% If no square root of preconditioner is available, sqrtPd = d (identity).
% Note that this may be incompatible with the preconditioner, if that one
% is supplied in the problem description. Always check with canGetPrecon
% and canGetSqrtPrecon.
%
% storedb is a StoreDB object, key is the StoreDB key to point x.
%
% See also: getPrecon canGetPrecon canGetSqrtPrecon getHessian

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, April 3, 2015.
% Contributors:
% Change log:

    % Allow omission of the key, and even of storedb.
    if ~exist('key', 'var')
        if ~exist('storedb', 'var')
            storedb = StoreDB();
        end
        key = storedb.getNewKey();
    end

    if isfield(problem, 'sqrtprecon')
    %% Apply sqrtprecon for the square root of the preconditioner

        % Check whether this function wants to deal with storedb or not.
        switch nargin(problem.sqrtprecon)
            case 2
                sqrtPd = problem.sqrtprecon(x, d);
            case 3
                % Obtain, pass along, and save the store for x.
                store = storedb.getWithShared(key);
                [sqrtPd, store] = problem.sqrtprecon(x, d, store);
                storedb.setWithShared(store, key);
            case 4
                % Pass along the whole storedb (by reference), with key.
                sqrtPd = problem.sqrtprecon(x, d, storedb, key);
            otherwise
                up = MException('manopt:getSqrtPrecon:badsqrtprecon', ...
                    'sqrtprecon should accept 2, 3 or 4 inputs.');
                throw(up);
        end

    else
    %% No preconditioner square root provided, so just use the identity.

        sqrtPd = d;

    end

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function out = MF_StateSpace_n4sid(y,ord,ptrain,steps)
% MF_StateSpace_n4sid   State space time-series model fitting.
%
% First fits the model to the whole time series, then trains it on the first
% portion and tries to predict the rest.
%
% In the second portion of this code, the state space model is fitted to the
% first p*N samples of the time series, where p is a given proportion and N is
% the length of the time series.
%
% This model is then used to predict the latter portion of the time
% series (i.e., the subsequent (1-p)*N samples).
%
% Model of the form:
% dx/dt = A x(t) + B u(t) + K e(t)
% y(t) = C x(t) + D u(t) + e(t)
% (for state space matrices A, B, C, D), disturbance matrix K (coefficients of
% noise input), input u, output y, vector of states x, and disturbance (noise) e.
%
%---INPUTS:
% y, the input time series
% ord, the order of state-space model to implement (can also be the string 'best')
% ptrain, the proportion of the time series to use for training
% steps, the number of steps ahead to predict
%
%---OUTPUTS: parameters from the model fitted to the entire time series, and
% goodness of fit and residual analysis from n4sid prediction.
%
% Uses the functions iddata, n4sid, aic, and predict from Matlab's System
% Identification Toolbox

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Check that a System Identification Toolbox license is available:
% ------------------------------------------------------------------------------
BF_CheckToolbox('identification_toolbox')

% ------------------------------------------------------------------------------
%% Check Inputs:
% ------------------------------------------------------------------------------
% (1) y: the time series as a column vector
% Convert y to time series object
N = length(y); % length of time series, N

y = iddata(y,[],1);

% (2) Order, the order of the state space model to fit. Can specify a positive
% integer or the string 'best'.
if nargin < 2 || isempty(ord)
    ord = 2;
end

% (3) train model on this proportion.
if nargin < 3 || isempty(ptrain)
    ptrain = 0.5;
end

% (4) steps, step-ahead prediction
if nargin < 4 || isempty(steps)
    steps = 1; % one-step ahead prediction
end

% ------------------------------------------------------------------------------
%% Build the state-space model
% ------------------------------------------------------------------------------
% Use the whole time series -- prediction comes later...
m = n4sid(y,ord); % fits a state-space model of given order

if strcmp(ord,'best')
    % also return the best order as an output statistic
    out.bestorder = length(m.k);
end

% ------------------------------------------------------------------------------
%% Model parameters
% ------------------------------------------------------------------------------
% Analysis of model

% Parameters
m_as = m.A; % 'transition' matrix in underlying ss model, x
m_ks = m.K; % coefficients for noise input in ss model
m_cs = m.C; % coefficients for measurement function, y
m_x0 = m.X0; % initial condition
m_np = length(m.ParameterVector); % number of parameters fitted

% Output model parameters
allm_as = m_as(:);
for i = 1:length(allm_as)
    out.(sprintf('A_%u',i)) = allm_as(i);
end
for i = 1:length(m_ks)
    out.(sprintf('k_%u',i)) = m_ks(i);
end
for i = 1:length(m_cs)
    out.(sprintf('c_%u',i)) = m_cs(i);
end
out.x0mod = sqrt(sum(m_x0.^2));
out.np = m_np; % the number of parameters, only a useful output if not specified

% Transition interval (this should always be 1 in this case, from how we've
% defined our time series, so not a useful output to record):
out.m_Ts = m.Ts;

% Goodness of fit outputs
% Since noisevar, lossfunction, and fpe so highly correlated, the default
% hctsa library only measures fpe.
out.m_noisevar = m.NoiseVariance; % a scalar number, basically the fpe
out.m_lossfn = m.EstimationInfo.LossFcn; % basically the fpe
out.m_fpe = m.EstimationInfo.FPE;
out.m_aic = aic(m);

% ------------------------------------------------------------------------------
%% Prediction
% ------------------------------------------------------------------------------

% Select first portion of data for estimation
% This could be any portion, actually... Maybe could look at robustness of
% model to different training sets...
ytrain = y(1:floor(ptrain*N));
% ytest = y;
ytest = y(floor(ptrain*N):end); % overlap

% Train the model on just this portion
% mp = armax(ytrain, orders);
try
    mp = n4sid(ytrain, ord);
catch emsg
    error('Couldn''t fit the model to this time series: %s',emsg.message)
    % out = NaN; return
end

%-------------------------------------------------------------------------------
% Step-ahead predictions
%-------------------------------------------------------------------------------
% steps = 2; % predicts this many steps ahead
% Maybe look at trends across different prediction horizons...
yp = predict(mp, ytest, steps, 'init', 'e'); % across whole ytest dataset

% plot the two:
% plot(y,yp);

mresiduals = ytest.y - yp.y;

%-------------------------------------------------------------------------------
% Statistics on residuals
%-------------------------------------------------------------------------------
residout = MF_ResidualAnalysis(mresiduals);

% Convert these to local outputs in quick loop
fields = fieldnames(residout);
for k = 1:length(fields);
    out.(fields{k}) = residout.(fields{k});
end

out.ac1diff = abs(CO_AutoCorr(y.y,1,'Fourier')) - abs(CO_AutoCorr(mresiduals,1,'Fourier'));

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [MB] = GB2MB(GB)
% Convert computery things from gigabytes to megabytes.
% Chad A. Greene 2012
MB = GB*1024 ;


text:
function [Jul1,Jul2]=GPS2TCG(Jul1,Jul2)
%GPS2TT  Convert from the timescale used by the Global Positioning System
%           (GPS)  given as a two-part Julian date to geocentric coordinate
%           time (TCG), represented as a two-part Julian date.
%
%INPUTS:    Jul1, Jul2  Two parts of a Julian date given in GPS time. The
%                       units of the date are days. The full date is the
%                       sum of both terms. The date is broken into two
%                       parts to provide more bits of precision. It does
%                       not matter how the date is split.
%
%OUTPUTS:   Jul1, Jul2  The time as a Julian date in TCG.
%
%This function just calls GPS2TAU and TAI2TCG.
%
%Many temporal coordinate systems standards are compared in [1].
%
%REFERENCES:
%[1] D. F. Crouse, "An Overview of Major Terrestrial, Celestial, and
%    Temporal Coordinate Systems for Target Tracking," Formal Report, Naval
%    Research Laboratory, no. NRL/FR/5344--16-10,279, 10 Aug. 2016, 173
%    pages.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

[Jul1,Jul2]=GPS2TAI(Jul1,Jul2);
[Jul1,Jul2]=TAI2TCG(Jul1,Jul2);
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:
classdef TestGFTTDetector
    %TestGFTTDetector

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','tsukuba_l.png');
        kfields = {'pt', 'size', 'angle', 'response', 'octave', 'class_id'};
    end

    methods (Static)
        function test_detect_img
            obj = cv.GFTTDetector('MaxFeatures',1000);
            assert(obj.MaxFeatures == 1000);
            typename = obj.typeid();

            img = imread(TestGFTTDetector.im);
            kpts = obj.detect(img);
            validateattributes(kpts, {'struct'}, {'vector'});
            assert(all(ismember(TestGFTTDetector.kfields, fieldnames(kpts))));
        end

        function test_detect_imgset
            img = imread(TestGFTTDetector.im);
            imgs = {img, img};
            obj = cv.GFTTDetector();

            kpts = obj.detect(imgs);
            validateattributes(kpts, {'cell'}, {'vector', 'numel',numel(imgs)});
            cellfun(@(kpt) validateattributes(kpt, {'struct'}, {'vector'}), kpts);
            cellfun(@(kpt) assert(all(ismember(TestGFTTDetector.kfields, fieldnames(kpt)))), kpts);
        end

        function test_detect_mask
            img = imread(TestGFTTDetector.im);
            mask = zeros(size(img), 'uint8');
            mask(:,1:end/2) = 255;  % only search left half of the image
            obj = cv.GFTTDetector();
            kpts = obj.detect(img, 'Mask',mask);
            validateattributes(kpts, {'struct'}, {'vector'});
            assert(all(ismember(TestGFTTDetector.kfields, fieldnames(kpts))));
            xy = cat(1, kpts.pt);
            assert(all(xy(:,1) <= ceil(size(img,2)/2)));
        end

        function test_error_1
            try
                cv.GFTTDetector('foobar');
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function main
close all
problemNames  = {'DenoisingEinstein';'AcceleratedDenoisingEinstein'};

for iproblem = 1:length(problemNames)
    pName = problemNames{iproblem};
    sLoader      = SettingsLoader(pName);
    settings     = sLoader.settings;
    denoisingProblem = DenoisingProblem(settings);
    denoisingProblem.solve();
    pD = denoisingProblem.plottingData;
    figure(1)
    hold on
    plot(pD.cost)
    figure(2)
    hold on
    plot(pD.dualGap)
end

end


text:
function pos=readsnx(sta_name,fname)

k=0; pos=zeros(1,3);
fid=fopen(fname);
while ~feof(fid)
    line = fgetl(fid);
    [tmp,line]=strtok(line);%#ok
    if strcmp(tmp,'+SOLUTION/APRIORI')
        while 1
            line = fgetl(fid);
            [tmp1,line]=strtok(line);%#ok
            [tmp2,line]=strtok(line);%#ok
            [tmp3,line]=strtok(line);%#ok
            if strcmp(tmp3,sta_name)
                [tmp,line]=strtok(line);%#ok
                [tmp,line]=strtok(line);%#ok
                [tmp,line]=strtok(line);%#ok
                [tmp,line]=strtok(line);%#ok
                [tmp,line]=strtok(line);%#ok
                [tmp,line]=strtok(line);%#ok
                pos(k+1)=str2num(tmp);%#ok
                k=k+1;
            end
            if (strcmp(tmp2,'STAZ') && strcmp(tmp3,sta_name)) || strcmp(tmp1,'-SOLUTION/APRIORI')
                break;
            end
        end
        break;
    end
end
fclose(fid);

return

repo: GINav


text:
function displayNeurons(obj, ind, C2, folder_nm)
%% view all components and delete components manually. it overlaps the
%   neuon contours with correlation image and shows the corresponding
%   spatial/temporal components
%% input:
%   ind: vector, indices of components to be displayed, no bigger than the maximum
%       number of neurons
%   C2:  K*T matrix, another temporal component to be displayed together
%       with the esitmated C. usually it is C without deconvolution.
%   folder_nm: string, the folder to output images neuron by neuron.

%% Author: Pengcheng Zhou, Carnegie Mellon University, 2016

if ~exist('ind', 'var') || isempty(ind)
    % display all neurons if ind is not specified.
    ind = 1:size(obj.A, 2);
end
if ~exist('C2', 'var'); C2=[]; end

if exist('folder_nm', 'var')&&(~isempty(folder_nm))
    % create a folder to save resulted images
    save_img = true;
    cur_cd = cd();
    if ~exist(folder_nm, 'dir'); mkdir(folder_nm);
    else
        fprintf('The folder has been created and old results will be overwritten. \n');
    end
    cd(folder_nm);
else
    save_img = false;
end
%% delete neurons with too few pixels
obj.delete(sum(obj.A>0, 1)<max(obj.options.min_pixel, 1));

ind_del = false(size(ind));     % indicator of deleting neurons
ctr = obj.estCenter();       % estimate neurons center
gSiz = obj.options.gSiz;        % maximum size of a neuron
Cn = obj.Cn;                    % correlation image
if isempty(Cn)
    fprintf('Please assign obj.Cn with correlation image!\n');
    return;
end
if isempty(obj.Coor) || (size(obj.A, 2)~=length(obj.Coor))   % contours of the neuron has not been calculated
    figure;
    obj.Coor = obj.get_contours();
end
Coor = obj.Coor;        % contours of all extracted neurons

% time
T = size(obj.C, 2);
t = 1:T;
if ~isnan(obj.Fs)
    t = t/obj.Fs;
    str_xlabel = 'Time (Sec.)';
else
    str_xlabel = 'Frame';
end

% start displaying neurons
figure('position', [100, 100, 1024, 512]);
m = 1;
while and(m>=1, m<=length(ind))
    %% contours + correlation image
    subplot(221); cla;
    obj.image(Cn, [0,1]); hold on; colormap winter;
    axis equal off tight;
    for k=1:m
        % plot contour
        tmp_con = Coor{ind(k)};
        cont_del = (sum(tmp_con<=1, 1)>0);
        tmp_con(:, cont_del) = [];
        if isempty(tmp_con)
            plot(ctr(m, 2), ctr(m, 2));
        else
            if and(k<m, ~ind_del(k))
                plot(tmp_con(1, 1:end), tmp_con(2, 1:end), 'color','k', 'linewidth', 2);
            elseif k==m
                plot(tmp_con(1, 1:end), tmp_con(2, 1:end), 'r', 'linewidth', 3);
            end
        end

        title(sprintf('Neuron %d', ind(m)));
    end

    %% zoomed-in view
    subplot(222);
    imagesc(reshape(obj.A(:, ind(m)), obj.options.d1, obj.options.d2));
    colormap(gca, jet);
    axis equal; axis off;
    x0 = ctr(ind(m), 2);
    y0 = ctr(ind(m), 1);
    xlim(x0+[-gSiz, gSiz]*2);
    ylim(y0+[-gSiz, gSiz]*2);

    %% temporal components
    subplot(2,2,3:4);cla;
    if ~isempty(C2)
        plot(t, C2(ind(m), :)*max(obj.A(:, ind(m))), 'linewidth', 2); hold on;
        plot(t, obj.C(ind(m), :)*max(obj.A(:, ind(m))), 'r');
    else

        plot(t, obj.C(ind(m), :)*max(obj.A(:, ind(m))));
    end
    xlabel(str_xlabel);

    %% save images
    if save_img
        saveas(gcf, sprintf('neuron_%d.png', ind(m)));
        m = m+1;
    else
        fprintf('Neuron %d, keep(k, default)/delete(d)/split(s)/delete all(da)/backward(b)/end(e):    ', ind(m));
        tmp_option = input('', 's');
        if tmp_option=='d'
            ind_del(m) = true;
            m = m+1;
        elseif strcmpi(tmp_option, 'b')
            m = m-1;
        elseif strcmpi(tmp_option, 'da')
            ind_del(m:end) = true;
            break;
        elseif strcmpi(tmp_option, 'k')
            ind_del(m) = false;
            m= m+1;
        elseif strcmpi(tmp_option, 's')
            try
                subplot(222);
                temp = imfreehand();
                tmp_ind = temp.createMask();
                tmpA = obj.A(:, ind(m));
                obj.A(:, end+1) = tmpA.*tmp_ind(:);
                obj.C(end+1, :) = obj.C(ind(m), :);
                obj.A(:, ind(m)) = tmpA.*(1-tmp_ind(:));
                obj.S(end+1, :) = obj.S(ind(m), :);
                obj.C_raw(end+1, :) = obj.C_raw(ind(m), :);
                obj.P.kernel_pars(end+1, :) = obj.P.kernel_pars(ind(m), :);
            catch
                sprintf('the neuron was not split\n');
            end
        elseif strcmpi(tmp_option, 't')
            try
                subplot(222);
                temp = imfreehand();
                tmp_ind = temp.createMask();
                obj.A(:, ind(m)) = obj.A(:, ind(m)).*tmp_ind(:);
            catch
                sprintf('the neuron was not trimmed\n');
            end
        elseif strcmpi(tmp_option, 'e')
            break;
        else
            m = m+1;
        end
    end
end
if save_img
    cd(cur_cd);
else
    obj.delete(ind(ind_del));
    obj.Coor = obj.get_contours(0.9);
end
figure; obj.viewContours(obj.Cn, 0.8, 0); close;

author: flatironinstitute

repo: CaImAn-MATLAB


text:
classdef ExistProbCalculatorX < BaseX
% ExistProbCalculatorX class
%
% Summary of ExistProbCalculatorX:
% This class can be used to compute (recursively) the existence probability
% of tracks.
%
% ExistProbCalculatorX Properties:
%   TrackList
%   AssocWeightsMatrix
%   ProbOfDetection
%   ProbOfDeath
%   ProbOfGating
%
% ExistProbCalculatorX Methods:
%   + ExistProbCalculatorX - Constructor method
%   + calculate - Calculate the existence probabilities
%
% (+) denotes public properties/methods
%
% March 2018 Lyudmil Vladimirov, University of Liverpool.

    properties
       TrackList
       AssocWeightsMatrix
       ProbOfDetection
       ProbOfDeath
       ProbOfGating
    end

    methods
        function this = ExistProbCalculatorX(varargin)
        % EXISTPROBCALCULATORX Constructor method
        %
        % Parameters
        % ----------
        % ProbOfDetection: scalar
        %   A scalar value in the range [0,1] which specifies the probability
        %   of detection for all targets
        % ProbOfDeath: scalar
        %   A scalar value in the range [0,1] which specifies the probability
        %   of death for all targets
        % ProbOfGating: scalar
        %   A scalar value in the range [0,1] which specifies the probability
        %   of gating for all targets
        %
        % Usage
        % -----
        % * ExistProbCalculatorX() returns a ExistProbCalculatorX object handle.
        % * ExistProbCalculatorX(ProbOfDetection,ProbOfDeath,ProbOfGating) returns
        %   a ExistProbCalculatorX( object handle configured with the provided
        %   probability of detection ProbOfDetection, probability of death ProbOfDeath
        %   and probability of gating ProbOfGating. This is particularly
        %   useful when all targets are characterised by the same parameters.
        % * ExistProbCalculatorX(config) can also be used, where config is a
        %   structure with fields config.ProbOfDetection, config.ProbOfDeath and
        %   config.ProbOfGating.
        % * PhdTrackInitiatorX(___,Name,Value) instantiates an object
        %   handle, configured with additional options specified by one or
        %   more Name,Value pair arguments.

            if(nargin==0)
                return;
            end

            % First check to see if a structure was received
            if(nargin==1)
                if(isstruct(varargin{1}))
                    config = varargin{1};
                    this.ProbOfDetection = config.ProbOfDetection;
                    this.ProbOfDeath = config.ProbOfDeath;
                    this.ProbOfGating = config.ProbOfGating;
                end
                return;
            elseif(nargin==3)
                this.ProbOfDetection = varargin{1};
                this.ProbOfDeath = varargin{2};
                this.ProbOfGating = varargin{3};
                return;
            end

            % Otherwise, fall back to input parser
            parser = inputParser;
            parser.KeepUnmatched = true;
            parser.parse(varargin{:});
            config = parser.Unmatched;
            this.ProbOfDetection = config.ProbOfDetection;
            this.ProbOfDeath = config.ProbOfDeath;
            this.ProbOfGating = config.ProbOfGating;
        end

        function argout = calculate(this,varargin)
        % CALCULATE Calculates the probability(ies) of existence
        %
        % Parameters
        % ----------
        % varargin{1}: (1 x NumTracks) cell array or TrackX object
        %   Either, a cell array of TrackX objects or a single TrackX
        %   object.
        % varargin{2}: (NumTracks x Nm+1) or (1 x Nm+1) matrix
        %   In the case that varargin{1} is a cell array, varargin{2}
        %   should be a (NumTracks x Nm+1) matrix, else it should be a
        %   (1 x Nm+1) row vector, where in both cases, each column contains
        %   the association weight between a measurement and a track. The
        %   first column corresponds to the dummy measurement.
        %
        % Returns
        % -------
        % argout: (1 x NumTracks) cell array or TrackX object
        %   Either, a cell array of TrackX objects or a single TrackX
        %   object, with their/its ProbOfExist field updated
        %
        % Usage
        % -----
        % * TrackList = calculate(this,TrackList,AssocWeightsMatrix)
        %   computes the new probability of existence for each TrackX
        %   object in TrackList and returns the updated TrackList
        % * Track = calculate(this,Track,AssocWeightsMatrix)
        %   computes the new probability of existence for the TrackX
        %   object Track and returns the updated Track
            if(isempty(varargin{1}))
                argout = varargin{1};
                return
            end
            if(iscell(varargin{1}))
                % Compute Existence Probabilities
                this.TrackList = varargin{1};
                this.AssocWeightsMatrix = varargin{2};
                numTracks = numel(this.TrackList);
                for trackInd = 1:numTracks
                    this.TrackList{trackInd}.ProbOfExist = (1 - this.ProbOfDeath)*this.TrackList{trackInd}.ProbOfExist;
                    denom = sum(this.AssocWeightsMatrix(trackInd,:))*this.TrackList{trackInd}.ProbOfExist + this.AssocWeightsMatrix(trackInd,1)*(1-this.TrackList{trackInd}.ProbOfExist)/((1-this.ProbOfDetection*this.ProbOfGating));
                    this.TrackList{trackInd}.ProbOfExist = (sum(this.AssocWeightsMatrix(trackInd,:))*this.TrackList{trackInd}.ProbOfExist)/denom;
                end
                argout = this.TrackList;
            elseif(isa(varargin{1},'TrackX'))
                Track = varargin{1};
                AssocWeights = varargin{2};
                Track.ProbOfExist = (1 - this.ProbOfDeath)*Track.ProbOfExist;
                denom = sum(AssocWeights)*Track.ProbOfExist + AssocWeights(1)*(1-Track.ProbOfExist)/((1-this.ProbOfDetection*this.ProbOfGating));
                Track.ProbOfExist = (sum(AssocWeights)*Track.ProbOfExist)/denom;
                argout = Track;
            end
        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [ r, s, area ] = node_reference_t3 ( )

%*****************************************************************************80
%
%% NODE_REFERENCE_T3 returns the basis nodes for the 3 node triangle.
%
%  Reference Element T3:
%
%    |
%    1  3
%    |  |\
%    |  | \
%    S  |  \
%    |  |   \
%    |  |    \
%    0  1-----2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real R(3), S(3), the coordinates of the basis nodes.
%
%    Output, real AREA, the area of the element.
%
  r(1:3) = [ 0.0, 1.0, 0.0 ];
  s(1:3) = [ 0.0, 0.0, 1.0 ];

  area = 0.5;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%--------------------------------------------------------------------------------------------------------
% The system is created based on the principles described in the following paper
% Jimmy SJ. Ren, Li Xu, Qiong Yan, Wenxiu Sun, "Shepard Convolutional Neural Networks",
% Advances in Neural Information Processing Systems (NIPS 2015)
% email: jimmy.sj.ren@gmail.com
%--------------------------------------------------------------------------------------------------------
addpath applications/Shepard_CNN/Shepard_super_res/
addpath applications/Shepard_CNN/utility/
addpath applications/image_denoise/utility/
addpath utils/
addpath cuda/
addpath mem/
addpath layers/
addpath layers_adapters/
addpath pipeline/

clearvars -global config;
clearvars -global mem;
clear gen_mask_patch_cat_idx_for_super_res;
clear;
global config;

% set5
% -- baby bird butterfly head woman
I = im2double(imread('applications/Shepard_CNN/Shepard_super_res/images/Set5/butterfly_GT.bmp'));
% set 14
% -- baboon barbara bridge coastguard comic face flowers foreman lenna man monarch pepper ppt3 zebra
I = im2double(imread('applications/Shepard_CNN/Shepard_super_res/images/Set14/ppt3.bmp'));
if(size(I,3) == 3)
    I = rgb2ycbcr(I);
end

I = modcrop(I, 4);
I = padarray(I, [4 4], 'replicate');

% downsample the image
D3chs = imresize(I, [size(I,1)/2 size(I,2)/2], 'bicubic');
D = D3chs(:,:,1);

% super-resolution x2
U = imresize(D, [size(I,1) size(I,2)], 'nearest');
prepare_sr_net(size(U, 1), size(U, 2), 'applications/Shepard_CNN/Shepard_super_res/results/shepard_layer_x2/sr_x2.mat');
mask = config.NEW_MEM([1 0; ...
                       0 0]);
mask = repmat(mask, size(U, 1)/2, size(U, 2)/2, 1);

final_output = apply_net_with_mask(config.NEW_MEM(U) .* mask, mask);
final_output = gather(final_output);
D3chs = imresize(D3chs, [size(final_output,1),size(final_output,2)]);
final_output = shave(final_output, [2 2]);
D3chs = shave(D3chs, [2 2]);

final_output = cat(3, final_output, D3chs(:,:,2:3));
final_output = ycbcr2rgb(double(final_output));
I_bicubic = ycbcr2rgb(D3chs);

figure('name', 'bicubic'), imshow(I_bicubic);
figure('name', 'Sheperd CNN'), imshow(final_output);
drawnow();

author: jimmy-ren