text:
%MSTRAJ Multi-segment multi-axis trajectory
%
% TRAJ = MSTRAJ(WP, QDMAX, TSEG, Q0, DT, TACC, OPTIONS) is a trajectory
% (KxN) for N axes moving simultaneously through M segment.  Each segment
% is linear motion and polynomial blends connect the segments.  The axes
% start at Q0 (1xN) and pass through M-1 via points defined by the rows of
% the matrix WP (MxN), and finish at the point defined by the last row of WP.
% The  trajectory matrix has one row per time step, and one column per
% axis.  The number of steps in the trajectory K is a function of the
% number of via points and the time or velocity limits that apply.
%
% - WP (MxN) is a matrix of via points, 1 row per via point, one column
%   per axis.  The last via point is the destination.
% - QDMAX (1xN) are axis speed limits which cannot be exceeded,
% - TSEG (1xM) are the durations for each of the K segments
% - Q0 (1xN) are the initial axis coordinates
% - DT is the time step
% - TACC (1x1) is the acceleration time used for all segment transitions
% - TACC (1xM) is the acceleration time per segment, TACC(i) is the acceleration
%   time for the transition from segment i to segment i+1.  TACC(1) is also
%   the acceleration time at the start of segment 1.
%
% TRAJ = MSTRAJ(WP, QDMAX, TSEG, [], DT, TACC, OPTIONS) as above but the
% initial coordinates are taken from the first row of WP.
%
% TRAJ = MSTRAJ(WP, QDMAX, Q0, DT, TACC, QD0, QDF, OPTIONS) as above
% but additionally specifies the initial and final axis velocities (1xN).
%
% Options::
% 'verbose'    Show details.
%
% Notes::
% - Only one of QDMAX or TSEG can be specified, the other is set to [].
% - If no output arguments are specified the trajectory is plotted.
% - The path length K is a function of the number of via points, Q0, DT
%   and TACC.
% - The final via point P(end,:) is the destination.
% - The motion has M segments from Q0 to P(1,:) to P(2,:) ... to P(end,:).
% - All axes reach their via points at the same time.
% - Can be used to create joint space trajectories where each axis is a joint
%   coordinate.
% - Can be used to create Cartesian trajectories where the "axes"
%   correspond to translation and orientation in RPY or Euler angle form.
% - If qdmax is a scalar then all axes are assumed to have the same
%   maximum speed.
%
% See also MTRAJ, LSPB, CTRAJ.

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

function [TG, t, info]  = mstraj(segments, qdmax, tsegment, q0, dt, Tacc, varargin)

    if isempty(q0)
        q0 = segments(1,:);
        segments = segments(2:end,:);
    end

    assert(size(segments,2) == size(q0,2), 'RTB:mstraj:badarg', 'WP and Q0 must have same number of columns');
    assert(xor(~isempty(qdmax), ~isempty(tsegment)), 'RTB:mstraj:badarg', 'Must specify either qdmax or tsegment, but not both');
    if isempty(qdmax)
        assert(length(tsegment) == size(segments,1), 'RTB:mstraj:badarg', 'Length of TSEG does not match number of segments');
    end
    if isempty(tsegment)
        if length(qdmax) == 1
            % if qdmax is a scalar assume all axes have the same speed
            qdmax = repmat(qdmax, 1, numcols(segments));
        end
        assert(length(qdmax) == size(segments,2), 'RTB:mstraj:badarg', 'Length of QDMAX does not match number of axes');
    end

    ns = numrows(segments);
    nj = numcols(segments);

    [opt,args] = tb_optparse([], varargin);

    if length(args) > 0
        qd0 = args{1};
    else
        qd0 = zeros(1, nj);
    end
    if length(args) > 1
        qdf = args{2};
    else
        qdf = zeros(1, nj);
    end

    % set the initial conditions
    q_prev = q0;
    qd_prev = qd0;

    clock = 0;      % keep track of time
    arrive = [];    % record planned time of arrival at via points

    tg = [];
    taxis = [];

    for seg=1:ns
        if opt.verbose
            fprintf('------------------- segment %d\n', seg);
        end

        % set the blend time, just half an interval for the first segment

        if length(Tacc) > 1
            tacc = Tacc(seg);
        else
            tacc = Tacc;
        end

        tacc = ceil(tacc/dt)*dt;
        tacc2 = ceil(tacc/2/dt) * dt;
        if seg == 1
            taccx = tacc2;
        else
            taccx = tacc;
        end

        % estimate travel time
        %    could better estimate distance travelled during the blend
        q_next = segments(seg,:);    % current target
        dq = q_next - q_prev;    % total distance to move this segment

        %% probably should iterate over the next section to get qb right...
        % while 1
        %   qd_next = (qnextnext - qnext)
        %   tb = abs(qd_next - qd) ./ qddmax;
        %   qb = f(tb, max acceleration)
        %   dq = q_next - q_prev - qb
        %   tl = abs(dq) ./ qdmax;

        if ~isempty(qdmax)
            % qdmax is specified, compute slowest axis

            qb = taccx * qdmax / 2;        % distance moved during blend
            tb = taccx;

            % convert to time
            tl = abs(dq) ./ qdmax;
            %tl = abs(dq - qb) ./ qdmax;
            tl = ceil(tl/dt) * dt;

            % find the total time and slowest axis
            tt = tb + tl;
            [tseg,slowest] = max(tt);

            info(seg).slowest = slowest;
            info(seg).segtime = tseg;
            info(seg).axtime = tt;
            info(seg).clock = clock;

            % best if there is some linear motion component
            if tseg <= 2*tacc
                tseg = 2 * tacc;
            end
        elseif ~isempty(tsegment)
            % segment time specified, use that
            tseg = tsegment(seg);
            slowest = NaN;
        end

        % log the planned arrival time
        arrive(seg) = clock + tseg;
        if seg > 1
            arrive(seg) = arrive(seg) + tacc2;
        end

        if opt.verbose
            fprintf('seg %d, slowest axis %d, time required %.4g\n', ...
                seg, slowest, tseg);
        end

        %% create the trajectories for this segment

        % linear velocity from qprev to qnext
        qd = dq / tseg;

        % add the blend polynomial
        qb = jtraj(q0, q_prev+tacc2*qd, 0:dt:taccx, qd_prev, qd);
        tg = [tg; qb(2:end,:)];

        clock = clock + taccx;     % update the clock

        % add the linear part, from tacc/2+dt to tseg-tacc/2
        for t=tacc2+dt:dt:tseg-tacc2
            s = t/tseg;
            q0 = (1-s) * q_prev + s * q_next;       % linear step
            tg = [tg; q0];
            clock = clock + dt;
        end

        q_prev = q_next;    % next target becomes previous target
        qd_prev = qd;
    end
    % add the final blend
    qb = jtraj(q0, q_next, 0:dt:tacc2, qd_prev, qdf);
    tg = [tg; qb(2:end,:)];
    info(seg+1).segtime = tacc2;
    info(seg+1).clock = clock;

    % plot a graph if no output argument
    if nargout == 0
        t = (0:numrows(tg)-1)'*dt;
        clf
        plot(t, tg, '-o');
        hold on
        plot(arrive, segments, 'bo', 'MarkerFaceColor', 'k');
        hold off
        grid
        xlabel('time');
        xaxis(t(1), t(end))
        return
    end
    if nargout > 0
        TG = tg;
    end
    if nargout > 1
        t = (0:numrows(tg)-1)'*dt;
    end
    if nargout > 2
        infout = info;
    end


text:
function dep = cfg_vout_gzip_files(job)

% Define virtual outputs for "Gzip Files". File names can either be
% assigned to a cfg_files input or to a evaluated cfg_entry.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_vout_gzip_files.m 6627 2015-12-04 08:20:49Z volkmar $

rev = '$Rev: 6627 $'; %#ok

dep            = cfg_dep;
dep.sname      = 'Gzipped Files';
dep.src_output = substruct('()',{':'});
dep.tgt_spec   = cfg_findspec({{'class','cfg_files', 'strtype','e'}});

author: spm


text:
function [C,phi,S12,S1,S2,f,zerosp,confC,phistd,Cerr]=coherencycpb(data1,data2,params,fscorr)
% Multi-taper coherency,cross-spectrum and individual spectra - continuous and binned point process data
%
% Usage:
%
% [C,phi,S12,S1,S2,f,zerosp,confC,phistd,Cerr]=coherencycpb(data1,data2,params,fscorr)
% Inputs:
%       data1 (continuous data in form samples x trials) -- required
%       data2 (binned spike data in form samples x trials) -- required
%       params: structure with fields tapers, pad, Fs, fpass, err, trialave
%       - optional
%           tapers : precalculated tapers from dpss or in the one of the following
%                    forms:
%                   (1) A numeric vector [TW K] where TW is the
%                       time-bandwidth product and K is the number of
%                       tapers to be used (less than or equal to
%                       2TW-1).
%                   (2) A numeric vector [W T p] where W is the
%                       bandwidth, T is the duration of the data and p
%                       is an integer such that 2TW-p tapers are used. In
%                       this form there is no default i.e. to specify
%                       the bandwidth, you have to specify T and p as
%                       well. Note that the units of W and T have to be
%                       consistent: if W is in Hz, T must be in seconds
%                       and vice versa. Note that these units must also
%                       be consistent with the units of params.Fs: W can
%                       be in Hz if and only if params.Fs is in Hz.
%                       The default is to use form 1 with TW=3 and K=5
%
%	        pad		    (padding factor for the FFT) - optional (can take values -1,0,1,2...).
%                    -1 corresponds to no padding, 0 corresponds to padding
%                    to the next highest power of 2 etc.
%			      	 e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
%			      	 to 512 points, if pad=1, we pad to 1024 points etc.
%			      	 Defaults to 0.
%           Fs   (sampling frequency) - optional. Default 1.
%           fpass    (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and Fs/2
%           err  (error calculation [1 p] - Theoretical error bars; [2 p] - Jackknife error bars
%                                   [0 p] or 0 - no error bars) - optional. Default 0.
%           trialave (average over trials when 1, don't average when 0) - optional. Default 0
%       fscorr   (finite size corrections, 0 (don't use finite size corrections) or 1
%              (use finite size corrections) - optional (available only for spikes). Defaults 0.
% Outputs:
%       C (magnitude of coherency - frequencies x trials if trialave=0; dimension frequencies if trialave=1)
%       phi (phase of coherency - frequencies x trials if trialave=0; dimension frequencies if trialave=1)
%       S12 (cross spectrum -  frequencies x trials if trialave=0; dimension frequencies if trialave=1)
%       S1 (spectrum 1 - frequencies x trials if trialave=0; dimension frequencies if trialave=1)
%       S2 (spectrum 2 - frequencies x trials if trialave=0; dimension frequencies if trialave=1)
%       f (frequencies)
%       zerosp (1 for trials where no spikes were found, 0 otherwise)
%       confC (confidence level for C at 1-p %) - only for err(1)>=1
%       phistd - theoretical/jackknife (depending on err(1)=1/err(1)=2) standard deviation for phi
%                Note that phi + 2 phistd and phi - 2 phistd will give 95% confidence
%                bands for phi - only for err(1)>=1
%       Cerr  (Jackknife error bars for C - use only for Jackknife -  err(1)=2)

if nargin < 2; error('Need data1 and data2'); end
if nargin < 3; params=[]; end;
[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
clear params
if nargin < 4 || isempty(fscorr); fscorr=0; end;

if nargout > 7 && err(1)==0;
%Cannot compute error bars if err(1)=0. Need to change params and run again.
    error('When errors are desired, err(1) has to be non-zero.');
end;
if nargout > 9 && err(1)~=2;
    error('Cerr computed only for Jackknife. Correct inputs and run again');
end;
[N,Ch]=check_consistency(data1,data2);
zerosp=zeros(1,Ch); % intialize the zerosp variable

nfft=max(2^(nextpow2(N)+pad),N);
[f,findx]=getfgrid(Fs,nfft,fpass);
tapers=dpsschk(tapers,N,Fs); % check tapers
J1=mtfftc(data1,tapers,nfft,Fs); % fourier transform of continuous data
[J2,Msp2,Nsp2]=mtfftpb(data2,tapers,nfft); % fourier transform of point process data
zerosp(Nsp2==0)=1; % set the trials where no spikes were found to have zerosp=1
J1=J1(findx,:,:);
J2=J2(findx,:,:);
S12=squeeze(mean(conj(J1).*J2,2));
S1=squeeze(mean(conj(J1).*J1,2));
S2=squeeze(mean(conj(J2).*J2,2));
if trialave; S12=squeeze(mean(S12,2)); S1=squeeze(mean(S1,2)); S2=squeeze(mean(S2,2)); end;
C12=S12./sqrt(S1.*S2);
C=abs(C12);
phi=angle(C12);
if nargout==10;
  if fscorr==1;
     [confC,phistd,Cerr]=coherr(C,J1,J2,err,trialave,[],Nsp2);
  else
     [confC,phistd,Cerr]=coherr(C,J1,J2,err,trialave);
  end;
elseif nargout==9;
  if fscorr==1;
     [confC,phistd]=coherr(C,J1,J2,err,trialave,[],Nsp2);
  else
     [confC,phistd]=coherr(C,J1,J2,err,trialave);
  end;
end;
clear Msp2

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [] = defaultlinecreatefcn(varargin)
% This is a createfcn for line objects.  To ensure that all lines created
% in the future automatically use LINECMENU, follow these instructions:
%
% 1.  Put this file and LINECMENU on your path.  I.e., put this file and
%     LINECMENU in one of the folders you see when you type:  path
% 2.  In STARTUP put this line:
%
%                 set(0,'defaultlinecreat',@defaultlinecreatefcn)
%
%    If you don't already have a STARTUP file, create one and make sure it
%    is on your path.  This is the file which MATLAB runs when it first
%    starts.
%
% The next time MATLAB is started, all lines will automatically have the
% LINECMENU options.  To make the changes effective without restarting
% MATLAB, simply execute the above call to SET from the command line.

set(varargin{1},'uicontextmenu',linecmenu);

author: Sable

repo: mcbench-benchmarks


text:
% jointprob() - rejection of odd columns of a data array  using
%              joint probability of the values in that column (and
%              using the probability distribution of all columns).
%
% Usage:
%   >>  [jp rej] = jointprob( signal );
%   >>  [jp rej] = jointprob( signal, threshold, jp, normalize, discret);
%
%
% Inputs:
%   signal     - one dimensional column vector of data values, two
%                dimensional column vector of values of size
%                sweeps x frames or three dimensional array of size
%                component x sweeps x frames. If three dimensional,
%                all components are treated independently.
%   threshold  - Absolute threshold. If normalization is used then the
%                threshold is expressed in standard deviation of the
%                mean. 0 means no threshold.
%   jp         - pre-computed joint probability (only perform thresholding).
%                Default is the empty array [].
%   normalize  - 0 = do not not normalize entropy. 1 = normalize entropy.
%                2 is 20% trimming (10% low and 10% high) proba. before
%                normalizing. Default is 0.
%   discret    - discretization variable for calculation of the
%                discrete probability density. Default is 1000 points.
%
% Outputs:
%   jp         - normalized joint probability  of the single trials
%                (size component x sweeps)
%   rej        - rejected matrix (0 and 1, size comp x sweeps)
%
% Remark:
%   The exact values of joint-probability depend on the size of a time
%   step and thus cannot be considered as absolute.
%
% See also: realproba()

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [jp, rej] = jointprob( signal, threshold, oldjp, normalize, discret );

if nargin < 1
	help jointprob;
	return;
end;
if nargin < 2
	threshold = 0;
end;
if nargin < 3
	oldjp = [];
end;
if nargin < 4
	normalize = 0;
end;
if nargin < 5
	discret = 1000;
end;

if size(signal,2) == 1 % transpose if necessary
	signal = signal';
end;

[nbchan pnts sweeps] = size(signal);
jp  = zeros(nbchan,sweeps);

if exist('oldjp') & ~isempty( oldjp ) % speed up the computation
	jp = oldjp;
else
	for rc = 1:nbchan

		% COMPUTE THE DENSITY FUNCTION
		% ----------------------------
		[ dataProba sortbox ] = realproba( signal(rc, :), discret );

		% compute all entropy
		% -------------------
		for index=1:sweeps
			datatmp = dataProba((index-1)*pnts+1:index*pnts);
			jp(rc, index) = - sum( log( datatmp ) );
			     % - sum( datatmp .* log( datatmp ) ); would be the entropy
		end;
	end;

	% normalize the last dimension
	% ----------------------------
	if normalize
        tmpjp = jp;
        if normalize == 2,
            tmpjp = sort(jp);
            tmpjp = tmpjp(round(length(tmpjp)*0.1):end-round(length(tmpjp)*0.1));
        end;
        try,
            switch ndims( signal )
             case 2,	jp = (jp-mean(tmpjp)) / std(tmpjp);
             case 3,	jp = (jp-mean(tmpjp,2)*ones(1,size(jp,2)))./ ...
                  (std(tmpjp,0,2)*ones(1,size(jp,2)));
            end;
        catch, error('Error while normalizing'); end;
	end;
end

% reject
% ------
if threshold ~= 0
    if length(threshold) > 1
    	rej = (threshold(1) > jp) | (jp > threshold(2));
    else
    	rej = abs(jp) > threshold;
    end;
else
	rej = zeros(size(jp));
end;

return;

author: PatternRecognition

repo: OpenBMI


text:
function out = spm_run_setlevel(job)
% out = spm_run_setlevel(job)
% test to see how likely it is that an SPM statistical image is a random field.
% based on:
%  Set-level threshold-free tests on the intrinsic volumes of SPMs.
%   Barnes GR, Ridgway GR, Flandin G, Woolrich M, Friston K. Neuroimage. 2013
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% $Id: spm_run_setlevel.m 5554 2013-06-13 09:13:56Z gareth $

%-Load SPM.mat file
%--------------------------------------------------------------------------
max_resel_dimension=3;  %% for now
LOWTHRESH=7; %% t value
threshlevels=[-LOWTHRESH:0.2:LOWTHRESH];
PLOTON=0;
ESTZEROLKC=0; %% don't try and estimate the 0th LKC (as this is given by the topology of the mask)

SPM = [];
load(job.spmmat{:});
out.spmmat = job.spmmat;

cindex=job.cindex; %% contrast of interest

%% First we need to redo the residual images - these have normally been cleaned up.

Ic=NaN; %% adjust for everything
disp('Writing new residual images');
Vres = spm_write_residuals(SPM,Ic);

%%% Now to make standarized residual images
nSres=size(Vres,2); %% number of subjects/observations
DIM=Vres.dim;
%-Initialise standardised residual images
%--------------------------------------------------------------------------

M=Vres.mat;
VResI(1:nSres) = deal(struct(...
    'fname',   [],...
    'dim',     DIM,...
    'dt',      [spm_type('float64') spm_platform('bigend')],...
    'mat',     M,...
    'pinfo',   [1 0 0]',...
    'descrip', 'spm_spm:StandardisedResiduals'));

for i = 1:nSres
    VResI(i).fname   = [sprintf('ResI_%04d', i) spm_file_ext];
    VResI(i).descrip = sprintf('spm_spm:ResI (%04d)', i);
end
VResI = spm_data_hdr_write(VResI);

ResMS=spm_read_vols(SPM.VResMS); %% residual mean square ( ResSS divided by tr(RV))

ResSS=ResMS.*SPM.xX.trRV; %%  residual sum of squares
%% now write out standardized residuals for each subject/observation
outvol            = NaN(size(ResSS));
cmask=find(isfinite(ResSS));

for j=1:nSres
    res=spm_read_vols(Vres(j));
    outvol(cmask) = res(cmask)./(sqrt(ResSS(cmask)/SPM.xX.erdf)); % or xX.erdf
    VResI(j) = spm_data_write(VResI(j), outvol);
end

maskV=SPM.VM;

mask=spm_read_vols(maskV);
maskind=intersect(cmask,find(mask));

datatrial1d=zeros(nSres,length(maskind));

%% disp read back in standardised residuals
for t=1:nSres,
    [sigvoldata,XYZ]=spm_read_vols(VResI(t));
    datatrial1d(t,:)=sigvoldata(maskind);
end % for t

%% read in statistical image
[Teststat,XYZ]=spm_read_vols(SPM.xCon(cindex).Vspm);
test_STAT=SPM.xCon(cindex).STAT;

ec_R0=zeros(nSres+1,1); %% topology
alleuler2_spm=zeros(nSres+1,length(threshlevels));

disp(sprintf('Getting ECs over %d residual images and %d thresholds...',nSres,length(threshlevels)));

for tr=1:nSres+1, % last extra trial is for Teststat
    vol_data=ones(size(sigvoldata)).*NaN;
    if tr<=nSres,
        vol_data(maskind)=datatrial1d(tr,:);
    else
        vol_data(maskind)=Teststat(maskind);
    end;

    binpoints=zeros(size(vol_data));
    binpoints(find(vol_data))=1;
    R0 = spm_resels_vol(binpoints,SPM.xVol.FWHM);
    ec_R0(tr)=R0(1);
    %% now get EC at each threshold and each observation
    for threshind=1:length(threshlevels),

        useind=find(vol_data>=threshlevels(threshind));
        binpoints=zeros(size(vol_data));
        binpoints(useind)=1;
        R0 = spm_resels_vol(binpoints,SPM.xVol.FWHM);
        ec_spm=R0(1);
        ECcount(tr,threshind)=ec_spm;
    end; % for threshind
end; % for tr

if max(ec_R0)~=min(ec_R0)
    error('subjects have different topology');
end;
R0_set=ec_R0(1);

%% GET EC DENSITIES FOR INDIVIDUAL TRIAL RESIDUALS (Z DIST) AND FOR TEST OF INTEREST (independent of data)
allpju_trial=[]; %% independent of data
allpju_test=[]; %% independent of data
df=[1 nSres-size(SPM.xX.X,2)];
trial_STAT='Z';
fitind=find(abs(threshlevels)<=LOWTHRESH);

for threshind=1:length(fitind),
    %%% NEED TO CHECK THIS ___
    [ECperresel_trial]=spm_ECdensity(trial_STAT,threshlevels(fitind(threshind)),df); %% ACTUALLY THIS IS LKC density
    [ECperresel_test]=spm_ECdensity(test_STAT,threshlevels(fitind(threshind)),df);

    for d=0:3, %% Turn density estimates from EC per resel to EC per LKC unit
        ECperLKC_trial(d+1,:)=ECperresel_trial(d+1,:)./power(4*log(2),d/2);
        ECperLKC_test(d+1,:)=ECperresel_test(d+1,:)./power(4*log(2),d/2);
    end;

    allpju_trial(threshind,:)=ECperLKC_trial(1:max_resel_dimension+1)';
    allpju_test(threshind,:)=ECperLKC_test(1:max_resel_dimension+1)';
end;

%%% GET RESEL ESTIMATES BASED ON SMOOTHNESS
[RPV]=spm_read_vols(SPM.xVol.VRpv); %% resel per voxel image
allR=SPM.xVol.R; %% resel counts

reselVec=SPM.xVol.R;
LKCresel=[];
for d=0:max_resel_dimension, %%% or alternatively LKC estimate from the reselts
    LKCresel(d+1)=reselVec(d+1)*power(4*log(2),d/2);
end;

for tbase=1:nSres+2,
    allpju=[];
    Y=[];
    if tbase<=nSres,
        epochtype=0; %% residual from trial
    else
        if tbase==nSres+1, %% test statistic
            epochtype=1;
        else
            epochtype=2; %% mean of ECs over trials
        end;
    end;

    switch epochtype,
        case 0, %% single residual trial
            Y=squeeze(ECcount(tbase,fitind))'; %% get LKC based on single residual image for these N trials
            allpju=allpju_trial; %% density for trial
        case 1 % single t stat
            Y=squeeze(ECcount(nSres+1,fitind))'; %% get LKC based on single Teststat image for these N trials
            allpju=allpju_test; %% density for test
        case 2, %% mean of all trials
            Y=squeeze(mean(ECcount(1:nSres,fitind),1))'; %% get LKC based on average EC over trials in iteration k
            %allaverageY(k,:)=Y;
            allpju=allpju_trial; %% density for trial
    end;

    dimension_test=max_resel_dimension;

    %% LKC based on average EC through basic regression
    if ~ESTZEROLKC, %% do not estimate 0th LKC
        LKC0=R0_set; %% TAKE THIS AS FIXED
        Ydash=Y-LKC0*allpju(:,1);
        useLKCind=2:dimension_test+1;
        LKC_est=pinv(allpju(:,useLKCind))*Ydash; %% get least squares estimate of extra LKC coeffs
        LKC=[LKC0; LKC_est];  %% put them back together with 0th order

    else %% estimate allLKCs
        useLKCind=1:dimension_test+1;
        LKC=pinv(allpju(:,useLKCind))*Ydash; %% get least squares estimate of all LKC coeffs
    end; % if

    allLKCregress(tbase,:)=LKC; %% LKC based on fit to average EC (/ or EC of Teststat)

end; % for tbase

gY = [squeeze(allLKCregress(nSres+1,useLKCind)); squeeze(allLKCregress(1:nSres,useLKCind))];
%then the design matrix and contrast could be e.g.
gX = [1 0; [zeros(nSres, 1) ones(nSres, 1)]];
gC = [1 -1]';
[CVA] = spm_cva(gY,gX,[],gC); %% do multivariate test

%% get empirical mean and sd of EC over subjects / observations.

meanECtest_regtrial=mean(squeeze(allLKCregress(1:nSres,:)))*allpju_test'; %% unweighted based on mean Euler
sdECtest_regtrial=std(squeeze(allLKCregress(1:nSres,:)))*allpju_test'; %% unweighted based on mean Euler

%% estimate what EC profile should be based on smoothness of image
meanECtest_resel=LKCresel*allpju_test';

%% PLOT RESULTS

Fgraph  = spm_figure('GetWin','Graphics'); spm_figure('Clear',Fgraph);

h=plot(threshlevels,squeeze(ECcount(nSres+1,:,:)),'r-',threshlevels,meanECtest_regtrial,'b:',threshlevels,meanECtest_resel,'go');
set(h,'LineWidth',3);
hold on;
errorbar(threshlevels,meanECtest_regtrial,sdECtest_regtrial,'.b');
set(gca,'Fontsize',18);
set(gcf,'color','w');
set(h,'LineWidth',3);
hold on;

xlabel('threshold');
ylabel('EC');
legend(sprintf('Observered EC for %s field',test_STAT),...
    sprintf('Random %s field based on regression',test_STAT),...
    sprintf('Random %s field based on smoothness',test_STAT));

title(sprintf('Probability that this is a random field  p<%3.4f ',CVA.p));

%%% FINISHED

%-Move to the directory where the SPM.mat file is
%--------------------------------------------------------------------------
original_dir = pwd;
cd(fileparts(job.spmmat{:}));
cd(original_dir);

fprintf('Done\n')
return


text:
% created: Zoya Bylinskii, Jan 2016

% Given a desired sigma blur value, this computes the cut off frequency
% required for the Gaussian low pass filter in the Fourier domain.
% This function runs Antonio's gaussian computation.

function [BF, gf] = run_antonioGaussian(img, sigma)

[sn, sm, c]=size(img);
n=max([sn sm]);

fc = n*sqrt(log(2)/(2*(pi^2)*(sigma^2)));

[BF, gf]=antonioGaussian(img, fc);

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
% READEGILOCS - look up locations for EGI EEG dataset.
%
% Usage:
%   >> EEG = readegilocs(EEG);
%   >> EEG = readegilocs(EEG, fileloc);
%
% Inputs:
%   EEG        - EEGLAB data structure
%
% Optional input:
%   fileloc    - EGI channel location file
%
% Outputs:
%   EEG        - EEGLAB data structure with channel location structure
%
% Author: Arnaud Delorme, SCCN/UCSD
%
% See also: POP_READEGI, READEGI

% Copyright (C) 12 Nov 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function EEG = readegilocs(EEG, fileloc);

if nargin < 1
    help readegilocs;
    return;
end

% importing channel locations
% ---------------------------
found = 1;
if nargin < 2
    switch EEG.nbchan
        case { 32 33 }, fileloc = 'GSN-HydroCel-32.sfp';
        case { 64 65 }, fileloc = 'GSN65v2_0.sfp';
        case { 128 129 }, fileloc = 'GSN129.sfp';
        case { 256 257 }, fileloc = 'GSN-HydroCel-257.sfp';
        otherwise, found = 0;
    end
end
if found
    fprintf('EGI channel location automatically detected %s ********* WARNING please check that this the proper file\n', fileloc);
    if EEG.nbchan == 64 || EEG.nbchan == 65 || EEG.nbchan == 256 || EEG.nbchan == 257
        fprintf( [ 'Warning: this function assumes you have a 64-channel system Version 2\n' ...
                   '         if this is not the case, update the channel location with the proper file' ]);
    end
    % remove the last channel for 33 channels

    peeglab = fileparts(which('eeglab.m'));
    fileloc = fullfile(peeglab, 'sample_locs', fileloc);
    locs = readlocs(fileloc);
    locs(1).type = 'FID';
    locs(2).type = 'FID';
    locs(3).type = 'FID';
    locs(end).type = 'REF';
    if EEG.nbchan == 256 || EEG.nbchan == 257
        if EEG.nbchan == 256
            chaninfo.nodatchans = locs([end]);
            locs([end]) = [];
        end
    elseif mod(EEG.nbchan,2) == 0,
        chaninfo.nodatchans = locs([1 2 3 end]);
        locs([1 2 3 end]) = [];
    else
        chaninfo.nodatchans = locs([1 2 3]);
        locs([1 2 3]) = [];
    end; % remove reference
    chaninfo.filename = fileloc;
    EEG.chanlocs   = locs;
    EEG.urchanlocs = locs;
    EEG.chaninfo   = chaninfo;
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function out=erf(x)

precision=x(1).precision;
out_rval=cell(size(x));
out_ival=cell(size(x));

for ii=1:numel(x)
 imag=false;
 [xrval,xival]=getVals(x,ii);
 if hasimag(xival), imag=true; end
 if imag
  error('argument of erf function for mp objects must be real');
 else
  out_rval{ii}=mpfr_erf(precision,xrval);
 end
end % for ii=1:max(ex,
out=class(struct('rval',out_rval,...
                 'ival',out_ival,...
                 'precision',precision),'mp');

author: opencobra

repo: cobratoolbox


text:
function imdb = cnn_imagenet_setup_data(varargin)
% CNN_IMAGENET_SETUP_DATA  Initialize ImageNet ILSVRC CLS-LOC challenge data
%    This function creates an IMDB structure pointing to a local copy
%    of the ILSVRC CLS-LOC data.
%
%    The ILSVRC data ships in several TAR archives that can be
%    obtained from the ILSVRC challenge website. You will need:
%
%    ILSVRC2012_img_train.tar
%    ILSVRC2012_img_val.tar
%    ILSVRC2012_img_test.tar
%    ILSVRC2012_devkit.tar
%
%    Note that images in the CLS-LOC challenge are the same for the 2012, 2013,
%    and 2014 edition of ILSVRC, but that the development kit
%    is different. However, all devkit versions should work.
%
%    In order to use the ILSVRC data with these scripts, please
%    unpack it as follows. Create a root folder <DATA>, by default
%
%    data/imagenet12
%
%    (note that this can be a simlink). Use the 'dataDir' option to
%    specify a different path.
%
%    Within this folder, create the following hierarchy:
%
%    <DATA>/images/train/ : content of ILSVRC2012_img_train.tar
%    <DATA>/images/val/ : content of ILSVRC2012_img_val.tar
%    <DATA>/images/test/ : content of ILSVRC2012_img_test.tar
%    <DATA>/ILSVRC2012_devkit : content of ILSVRC2012_devkit.tar
%
%    In order to speedup training and testing, it may be a good idea
%    to preprocess the images to have a fixed size (e.g. 256 pixels
%    high) and/or to store the images in RAM disk (provided that
%    sufficient RAM is available). Reading images off disk with a
%    sufficient speed is crucial for fast training.

opts.dataDir = fullfile('data','imagenet12') ;
opts.lite = false ;
opts = vl_argparse(opts, varargin) ;

% -------------------------------------------------------------------------
%                                                  Load categories metadata
% -------------------------------------------------------------------------

d = dir(fullfile(opts.dataDir, '*devkit*')) ;
if numel(d) == 0
  error('Make sure that ILSVRC data is correctly installed in %s', ...
    opts.dataDir) ;
end
devkitPath = fullfile(opts.dataDir, d(1).name) ;

% find metadata
mt = dir(fullfile(devkitPath, 'data', 'meta_clsloc.mat')) ;
if numel(mt) == 0
  mt = dir(fullfile(devkitPath, 'data', 'meta.mat')) ;
end
metaPath = fullfile(devkitPath, 'data', mt(1).name) ;

% find validation images labels
tmp = dir(fullfile(devkitPath, 'data', '*_validation_ground_truth*')) ;
valLabelsPath = fullfile(devkitPath, 'data', tmp(1).name) ;

% find validation images blacklist
tmp = dir(fullfile(devkitPath, 'data', '*_validation_blacklist*')) ;
if numel(tmp) > 0
  valBlacklistPath = fullfile(devkitPath, 'data', tmp(1).name) ;
else
  valBlacklistPath = [] ;
  warning('Could not find validation images blacklist file');
end

fprintf('using devkit %s\n', devkitPath) ;
fprintf('using metadata %s\n', metaPath) ;
fprintf('using validation labels %s\n', valLabelsPath) ;
fprintf('using validation blacklist %s\n', valBlacklistPath) ;

meta = load(metaPath) ;
cats = {meta.synsets(1:1000).WNID} ;
descrs = {meta.synsets(1:1000).words} ;

imdb.classes.name = cats ;
imdb.classes.description = descrs ;
imdb.imageDir = fullfile(opts.dataDir, 'images') ;

% -------------------------------------------------------------------------
%                                                           Training images
% -------------------------------------------------------------------------

fprintf('searching training images ...\n') ;
names = {} ;
labels = {} ;
for d = dir(fullfile(opts.dataDir, 'images', 'train', 'n*'))'
  [~,lab] = ismember(d.name, cats) ;
  ims = dir(fullfile(opts.dataDir, 'images', 'train', d.name, '*.JPEG')) ;
  names{end+1} = strcat([d.name, filesep], {ims.name}) ;
  labels{end+1} = ones(1, numel(ims)) * lab ;
  fprintf('.') ;
  if mod(numel(names), 50) == 0, fprintf('\n') ; end
  %fprintf('found %s with %d images\n', d.name, numel(ims)) ;
end
names = horzcat(names{:}) ;
labels = horzcat(labels{:}) ;

if numel(names) ~= 1281167
  warning('Found %d training images instead of 1,281,167. Dropping training set.', numel(names)) ;
  names = {} ;
  labels =[] ;
end

names = strcat(['train' filesep], names) ;

imdb.images.id = 1:numel(names) ;
imdb.images.name = names ;
imdb.images.set = ones(1, numel(names)) ;
imdb.images.label = labels ;

% -------------------------------------------------------------------------
%                                                         Validation images
% -------------------------------------------------------------------------

ims = dir(fullfile(opts.dataDir, 'images', 'val', '*.JPEG')) ;
names = sort({ims.name}) ;
labels = textread(valLabelsPath, '%d') ;

if numel(ims) ~= 50e3
  warning('Found %d instead of 50,000 validation images. Dropping validation set.', numel(ims))
  names = {} ;
  labels =[] ;
else
  if ~isempty(valBlacklistPath)
    black = textread(valBlacklistPath, '%d') ;
    fprintf('blacklisting %d validation images\n', numel(black)) ;
    keep = setdiff(1:numel(names), black) ;
    names = names(keep) ;
    labels = labels(keep) ;
  end
end

names = strcat(['val' filesep], names) ;

imdb.images.id = horzcat(imdb.images.id, (1:numel(names)) + 1e7 - 1) ;
imdb.images.name = horzcat(imdb.images.name, names) ;
imdb.images.set = horzcat(imdb.images.set, 2*ones(1,numel(names))) ;
imdb.images.label = horzcat(imdb.images.label, labels') ;

% -------------------------------------------------------------------------
%                                                               Test images
% -------------------------------------------------------------------------

ims = dir(fullfile(opts.dataDir, 'images', 'test', '*.JPEG')) ;
names = sort({ims.name}) ;
labels = zeros(1, numel(names)) ;

if numel(labels) ~= 100e3
  warning('Found %d instead of 100,000 test images', numel(labels))
end

names = strcat(['test' filesep], names) ;

imdb.images.id = horzcat(imdb.images.id, (1:numel(names)) + 2e7 - 1) ;
imdb.images.name = horzcat(imdb.images.name, names) ;
imdb.images.set = horzcat(imdb.images.set, 3*ones(1,numel(names))) ;
imdb.images.label = horzcat(imdb.images.label, labels) ;

% -------------------------------------------------------------------------
%                                                            Postprocessing
% -------------------------------------------------------------------------

% sort categories by WNID (to be compatible with other implementations)
[imdb.classes.name,perm] = sort(imdb.classes.name) ;
imdb.classes.description = imdb.classes.description(perm) ;
relabel(perm) = 1:numel(imdb.classes.name) ;
ok = imdb.images.label >  0 ;
imdb.images.label(ok) = relabel(imdb.images.label(ok)) ;

if opts.lite
  % pick a small number of images for the first 10 classes
  % this cannot be done for test as we do not have test labels
  clear keep ;
  for i=1:10
    sel = find(imdb.images.label == i) ;
    train = sel(imdb.images.set(sel) == 1) ;
    val = sel(imdb.images.set(sel) == 2) ;
    train = train(1:256) ;
    val = val(1:40) ;
    keep{i} = [train val] ;
  end
  test = find(imdb.images.set == 3) ;
  keep = sort(cat(2, keep{:}, test(1:1000))) ;
  imdb.images.id = imdb.images.id(keep) ;
  imdb.images.name = imdb.images.name(keep) ;
  imdb.images.set = imdb.images.set(keep) ;
  imdb.images.label = imdb.images.label(keep) ;
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
 function h = sgtitlef(varargin)
%function h = sgtitlef(varargin)
%|
%| version of sgtitle with built-in sprintf
%| and defaults to latex interpreter
%| also supports default font size from ir_fontsize()
%|
%| Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end

opt = {'fontsize', ir_fontsize('title')};
%opt = {opt{:}, 'fontname', 'times'};

if ir_is_octave
	opt = {opt{:}, 'fontname', 'Helvetica'};
end

varargin{1} = ir_strrep_tex(varargin{1}); % trick

if ir_is_octave
	tex = {};
	str = varargin{1};
	str = strrep(str, '$', ''); % remove $ needed for latex
	varargin{1} = str;
else
	tex = {'interpreter', 'latex'};
	str = varargin{1};
	str = strrep(str, '\', '\\'); % because of sprintf below
	% trick to replace "|" with "$|$" in strings
	if ~isempty(strfind(str, '|')) && isempty(strfind(str, '$'))
		str = strrep(str, '|', '$|$');
	end
	varargin{1} = str;
end

if isfreemat
	for ii=1:length(varargin)
		if streq(varargin{ii}, 'interpreter') % not supported by freemat
			varargin{ii} = {};
			varargin{ii+1} = {};
		end
	end
end

if im
	tmp = sprintf(varargin{:});
	tmp = strrep(tmp, '%', '\%'); % else latex will ignore %
	hh = sgtitle(tmp, tex{:}, opt{:});
else
	hh = [];
end

if nargout
	h = hh;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% Compute the center frequency of the mel filterbank
% Inputs:
%   start_freq:     the lowerest linear frequency included in the calculation
%   linear_samp:    sampling frequency
%   N_mel:          number of mel banks used
% Output:
%   mel_center_freq:    the center linear frequency of each mel bank
%   bin_upper_edge:

function [mel_center_freq, bin_upper_edge] = mel_center_FE(start_freq, linear_samp, N_mel)

% get the mel version of start frequency and linear sampling rate
start_mel = linear2mel(start_freq);
samp_mel = linear2mel(linear_samp/2);

for i=1:N_mel
    tmp_mel = (i-1)*(samp_mel-start_mel) / (N_mel+1);
    mel_center_freq(i) = mel2linear(tmp_mel+start_mel);
    tmp_mel = (i+1)*(samp_mel-start_mel) / (N_mel+1);
    bin_upper_edge(i) =  mel2linear(tmp_mel+start_mel);
end

author: singaxiong

repo: SignalGraph


text:
function outMap=computeCorrelationMap(view,scan,timeSeries,normalizeFlag)
% outMap=computeCorrelationMap(view,scan,timeSeries,normalizeFlag)
% Computes a parameter map that is the cross correlation of each voxel (in
% the inplance view - for the current dt and specified scan with the reference
% timeSeries.
%
% The procedure happens slice-by-slice and (normalized) cross correlation
% is computed in the Fourier domain to speed things up.
% The thing returned is the COMPLEX map of the zero-th lag.
% You can do abs(map)
% ARW 103106: Wrote it.
%

% Check the inputs
if (ieNotDefined('view'))
    error('You must supply an INPLANE view');
end
if (ieNotDefined('scan'))
    error('You must supply a scan index');
end
if (ieNotDefined('timeSeries'))
    error('You must supply a reference time series');
end
if (ieNotDefined('normalizeFlag'))
    warning('Doing a normalized cross correlation by default');
    normalizeFlag=1;
end

% Pre-allocate the output
dataSize=viewget(view,'datasize');
outMap=zeros(dataSize(1),dataSize(2),dataSize(3));

nSlices=dataSize(3);
nTRs=viewget(view,'nframes');

% Check the reference tSeries
timeSeries=timeSeries(:);
if length(timeSeries)~=nTRs
    error ('The reference time series is a different length to the data');
end
disp(scan);

% We're going to do the xcorr in fourier space
% That means that we need to pre-compute the fft of the ref tSeries
% and then make a matrix..

% Always remove the mean
timeSeries=timeSeries-mean(timeSeries(:));

fRefTSeries=fft(timeSeries(:));
fRefTSeries=conj(fRefTSeries(:));  %Compute the complex conjugate here so that we get xcorr in Fourier space

if (normalizeFlag) % Ensure max corr is 1
    fRefTSeries=fRefTSeries./sqrt(sum(abs(fRefTSeries(:)).^2));
end

nVoxPerSlice=dataSize(1)*dataSize(2); % Compute the number of voxels in the slice
FFT_tSeriesMatrix=repmat(fRefTSeries,1,nVoxPerSlice); % Make an array of identical fourier transforms of the input to speed multiplication. Freq goes down

for thisSlice=1:nSlices % loop over slices

    % Get the tSeries data for this slice
    thisTSeries=loadtSeries(view, scan,thisSlice);

    % Remove the mean
    thisTSeries=thisTSeries-repmat(mean(thisTSeries),nTRs,1);

    % We can do the normXCorr right away
    FFT_thisTSeries=fft(thisTSeries);

    if (normalizeFlag) % Normalize the power. Check this - there might be an additional factor of nframes in the final result
        FFT_thisTSeries=FFT_thisTSeries./repmat(sqrt(sum(abs(FFT_thisTSeries).^2)),nTRs,1);
    end

    xcorrTSeries=FFT_thisTSeries.*(FFT_tSeriesMatrix); % Complex conjugate of the reference already computed

    i_xcorr=ifft(xcorrTSeries); % Do the ifft to get the cross correlation
    compData=i_xcorr(1,:); % Take the zero-th lag. This assumes that there are no lagged copies of the signal.

    % Re-shape so that you can place it directly into a parameter map.
    outMap(:,:,thisSlice)=reshape(compData,dataSize(1),dataSize(2));

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Copyright (C) Daphne Koller, Stanford Univerity, 2012
%
% Author: Huayan Wang, Andrew Duchi

% visualize a configuration the body parts

function img = ShowPose( pose )

% pose 10 x 3. 10 : body parts,
%               3 : y, x, alpha

pose(:,1) = pose(:,1) + 100;
pose(:,2) = pose(:,2) + 150;

pose = reshape(pose, [10 3]);
part_length = [60, 20, 32, 33, 32, 33, 46, 49, 46, 49];
part_width = [18, 10, 7, 5, 7, 5, 10, 7, 10, 7];
img = zeros(300, 300);
for part = 1:10

    startpt = round(pose(part, 1:2));
    axis = [sin(pose(part,3) - pi/2) cos(pose(part,3) - pi/2)];
    xaxis = [cos(pose(part,3) - pi/2) -sin(pose(part,3) - pi/2)];
    endpt = round(startpt + part_length(part) * axis);

    corner1 = round(startpt + xaxis * part_width(part));
    corner2 = round(startpt - xaxis * part_width(part));
    corner3 = round(endpt + xaxis * part_width(part));
    corner4 = round(endpt - xaxis * part_width(part));

   img = func_DrawLine(img, corner1(1), corner1(2), corner2(1), corner2(2),1);
   img = func_DrawLine(img, corner1(1), corner1(2), corner3(1), corner3(2),1);
   img = func_DrawLine(img, corner4(1), corner4(2), corner2(1), corner2(2),1);
   img = func_DrawLine(img, corner4(1), corner4(2), corner3(1), corner3(2),1);

   if startpt(1) > 3 && startpt(1) < 298 && startpt(2) > 3 && startpt(2) < 298
       img(startpt(1)-3 : startpt(1) + 3, startpt(2) - 3 : startpt(2)+3) = ones(7,7);
   end
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [ kids, ntypelist, layers ] = roomhierarchy( data, RELATIONS, NODETYPES )
%ROOMHIERARCHY

obblist = data.obblist;
Rmat = data.Rmat;
Dmat = data.Dmat;
labellist = data.labellist;
sur = data.sur;
objnum = size(obblist,2);

%% build support tree
supportparent = zeros(size(obblist,2),1);
for i = 1:objnum
    rlist = Rmat(:,i);
    pp = find(rlist==RELATIONS.SUPPORT1);
    if(length(pp)>1)
        pd = Dmat(i,pp);
        [~,I] = min(pd);
        supportparent(i) = pp(I);
    end
    if(length(pp)==1)
        supportparent(i) = pp;
    end
end

%% compute the layer values
% nobjlist = length(supportparent);% add the floor
% layers = zeros(size(obblist,2),1);
% layers(nobjlist(1)) = 1;
% while(length(nobjlist)>0)
%     i = nobjlist(1);
%     nobjlist = nobjlist(2:end);
%     ind = find(supportparent==i);
%     nobjlist = [nobjlist;ind];
%     layers(ind) = layers(i)+1;
% end

%% split by surround relation
parent = [supportparent(:)';zeros(1,length(supportparent))+NODETYPES.SUPPORT];
flooridx = strmatch('floor',labellist,'exact');
parent(1,flooridx) = 0;
% a central object may be involved with multiple surround relations,
% but it can have only one parent
for i = 1:objnum
    % find the possible parents, sur pair id, distance in the relation
    pdlist = [];
    for j = 1:length(sur)
        slist = sur{j};
        if(slist(1)==i&&length(slist)>2)
            t = tabulate(supportparent(slist(2:end)));
            if(size(t,2)>1)
                [~,pp] = max(t(:,2));
            else
                pp = 1;
            end
            pd = [t(pp,1),j,Dmat(slist(1),slist(2))];
            pdlist = [pdlist;pd];
        end
    end
    surind = [];
    if(size(pdlist,1)>0)
        ind = find(pdlist(:,1)==parent(1,i));
        surind = pdlist(ind,2);
    end

    if(length(surind)>1)
        slist = sur{surind(1)};
        surchildlist = slist(2:end);
        innodelist = [];
        for j = 1:length(surchildlist)
            plist = [surchildlist(j)];
            for k = 2:length(surind)
                sklist = sur{surind(k)};
                sklist = sklist(2:end);
                dlist = Dmat(surchildlist(j),sklist);
                [~,I] = min(dlist);
                plist = [plist,sklist(I)];
            end
            parent(:,size(parent,2)+1) = [0;NODETYPES.COOCCUR]; % note here may make mistake
            parent(1,plist) = size(parent,2);
            obblist(:,size(obblist,2)+1) = updateOBB(obblist,plist);
            innodelist = [innodelist;size(parent,2)];
        end
%         parent(innodelist) = i;% assign the central obj as their parent
        innodelist = [i;innodelist(:)];
        parent(:,size(parent,2)+1) = [parent(1,i);NODETYPES.SURROUND];
        parent(1,innodelist) = size(parent,2); % assign the central obj as their parent
        obblist(:,size(obblist,2)+1) = updateOBB(obblist,innodelist);
    end

    if(length(surind)==1)
        slist = sur{surind(1)};
        parent(:,size(parent,2)+1) = [parent(1,i);NODETYPES.SURROUND];
        parent(1,slist) = size(parent,2);
        obblist(:,size(obblist,2)+1) = updateOBB(obblist,slist);
    end
end

%% create the distmat for all the current nodes
distmat = Dmat;
for i = 1+objnum:size(parent,2)
    ind = find(parent(1,:)==i);
    dm = distmat(ind,:);
    if(length(ind)>1)
        dm = min(dm);
    end
    if(length(ind)>0)
        distmat = [distmat; dm];
        dm = [dm(:);0];
        distmat = [distmat,dm];
    end
end

wallidx = strmatch('wall',labellist,'exact');
flooridx = strmatch('floor',labellist,'exact');
if(length(wallidx)>0)
    %% create wall nodes
    %% shoud be changed to obj cluster algorithm
    floorchildlist = find(parent(1,:)==flooridx);
    floorchildlist = setdiff(floorchildlist,wallidx);
    [ clusters ] = objclustering( obblist, distmat, wallidx, floorchildlist );
     for j = 1:length(clusters)
        list = clusters{j};
        parent(1,list) = wallidx(j);
    end
%     for i = 1:length(floorchildlist)
%         cid = floorchildlist(i);
%         distlist = distmat(cid,wallidx);
%         [~,I] = min(distlist);
%         parent(:,cid) = [wallidx(I);parent(2,cid)];
%     end
%
end
parent(2,wallidx) = NODETYPES.WALL;
parent(2,flooridx) = NODETYPES.ROOM;%%%%%%%%%

% colorlist = {'r','g','b','k'};
% for i = 1:length(clusters)
%     c = clusters{i};
%     draw3dOBB_v2(obblist(:,wallidx(i)),colorlist{i});
%     for j = 1:length(c)
%         draw3dOBB_v2(obblist(:,c(j)),colorlist{i});
%     end
% end

[ kids, ntypelist ] = organize_binarytree( parent, distmat, objnum, NODETYPES, flooridx);

%% check and correct the ntypelist
for i = 1:length(ntypelist)
    flag = ntypelist(i);
    if(flag==NODETYPES.COOCCUR)
        k = kids{i};
        if(k(1)<size(Rmat,1)&k(2)<size(Rmat,1))
            rel = Rmat(k(1),k(2));
            if(rel==RELATIONS.SUPPORT1)
                ntypelist(i) = NODETYPES.SUPPORT;
            end
            if(rel==RELATIONS.SUPPORT2)
                ntypelist(i) = NODETYPES.SUPPORT;
                tmp = k(1);
                k(1) = k(2);
                k(2) = tmp;
            end
            kids{i} = k;
        end
    end
end

end

repo: GRAINS


text:
      satis_distort = 0;

      while ~satis_distort,

          k_g_save = k_g;

	 k_g = input(['Guess for distortion factor kc ([]=' num2str(k_g_save) '): ']);

	 if isempty(k_g), k_g = k_g_save; end;

x_n = (x - c_g(1))/f_g;
y_n = (y - c_g(2))/f_g;

[x_pn] = comp_fisheye_distortion([x_n' ; y_n'],[k_g;0;0;0]);

% Compute the inside points through computation of the planar homography (collineation)

a00 = [x_pn(1,1);x_pn(2,1);1];
a10 = [x_pn(1,2);x_pn(2,2);1];
a11 = [x_pn(1,3);x_pn(2,3);1];
a01 = [x_pn(1,4);x_pn(2,4);1];

% Compute the planar collineation: (return the normalization matrix as well)
[Homo,Hnorm,inv_Hnorm] = compute_homography([a00 a10 a11 a01],[0 1 1 0;0 0 1 1;1 1 1 1]);

% Build the grid using the planar collineation:

x_l = ((0:n_sq_x)'*ones(1,n_sq_y+1))/n_sq_x;
y_l = (ones(n_sq_x+1,1)*(0:n_sq_y))/n_sq_y;
pts = [x_l(:) y_l(:) ones((n_sq_x+1)*(n_sq_y+1),1)]';

XXpn = Homo*pts;
XXpn = XXpn(1:2,:) ./ (ones(2,1)*XXpn(3,:));

XX = apply_fisheye_distortion(XXpn,[k_g;0;0;0]);

XX(1,:) = f_g*XX(1,:) + c_g(1);
XX(2,:) = f_g*XX(2,:) + c_g(2);

	 figure(2);
	 image(I);
	 colormap(map);
	 zoom on;
	 hold on;
	 plot(XX(1,:),XX(2,:),'r+');
	 title('The red crosses should be on the grid corners...');
	 hold off;

	 satis_distort = input('Satisfied with distortion? ([]=no, other=yes) ');

	 satis_distort = ~isempty(satis_distort);

      end;

repo: ekfmonoslam


text:
function test_ft_appendspike

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_appendspike

spike = [];
spike.label = {'unit1'  'unit2'  'unit3'};
spike.timestamp = {linspace(0,100,100) linspace(0,200,200) linspace(0,300,300)};
spike.trial = {randn(1,100) randn(1,200) randn(1,300)};
spike.trialtime = [linspace(0,30,300)', linspace(0,30,300)'+0.5];

spike2 = spike;
spike.label = {'unit12'  'unit22'  'unit32'};
spike2.trial = {randn(1,100) randn(1,200) randn(1,300)};

cfg = [];
cfg.trl = spike.trialtime;
spikeout = ft_appendspike(cfg,spike,spike2);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [maxVal, maxPos] = max(f)
%MAX   Global maximum of a CHEBTECH on [-1,1].
%   MAXVAL = MAX(F) returns the global maximum of the CHEBTECH F on [-1,1].  If
%   F is an array-valued CHEBTECH, MAXVAL is a row vector whose Kth entry is the
%   global maximum of the Kth column of F.
%
%   [MAXVAL, MAXPOS] = MAX(F) returns also a value such that MAXVAL = F(MAXPOS).
%
%   If F is complex-valued then absolute values are taken to determine maxima
%   but the resulting value corresponds to that of the original function. That
%   is, MAXVAL = feval(F, MAXPOS) where [~, MAXPOS] = MAX(abs(F));
%
% See also MIN, MINANDMAX.

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

% To avoid code duplication, we simply call MINANDMAX:
[maxVal, maxPos] = minandmax(f);
maxVal = maxVal(2,:);
maxPos = maxPos(2,:);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% EX_STOKES_SQUARE_SG: solve the Stokes problem in the unit square with the SubGrid method.

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

% Type of boundary conditions for each side of the domain
problem_data.drchlt_sides = 1:4;
problem_data.nmnn_sides = [];

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

% Force term
fx = @(x, y) (6 * x + y .* cos(x .* y) + 2 * cos(y) .* sin(x));
fy = @(x, y) (x .* cos(x .* y) - 2 * cos(x) .* sin(y));
problem_data.f  = @(x, y) cat(1, ...
                reshape (fx (x,y), [1, size(x)]), ...
                reshape (fy (x,y), [1, size(x)]));

% Boundary terms
uxex = @(x, y) (sin(x) .* cos(y));
uyex = @(x, y) (-sin(y) .* cos(x));

problem_data.h = @(x, y, iside) cat(1, ...
                  reshape (uxex (x,y), [1, size(x)]), ...
                  reshape (uyex (x,y), [1, size(x)]));

% Exact solution, to compute the errors
problem_data.velex = @(x, y) cat(1, ...
                  reshape (uxex (x,y), [1, size(x)]), ...
                  reshape (uyex (x,y), [1, size(x)]));

problem_data.pressex = @(x, y) (3 * x.^2 + sin(x .* y) - 1.239811742000564725943866);

problem_data.gradvelex = @test_stokes_square_bc_graduex;

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.element_name = 'SG';   % Element type for discretization
method_data.degree       = [ 3  3];  % Degree of the splines (pressure space)
method_data.regularity   = [ 2  2];  % Regularity of the splines (pressure space)
method_data.nsub         = [10 10];  % Number of subdivisions
method_data.nquad        = [ 5  5];  % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space_v, vel, space_p, press] = ...
                       solve_stokes (problem_data, method_data);

% 4) POST-PROCESSING
% 4.1) COMPARISON WITH EXACT SOLUTION
error_l2_p = sp_l2_error (space_p, msh, press, problem_data.pressex)
[error_h1_v, error_l2_v] = ...
   sp_h1_error (space_v, msh, vel, problem_data.velex, problem_data.gradvelex)

% 4.2) EXPORT TO PARAVIEW
output_file = 'SQUARE_SG_Deg3_Reg2_Sub10';

fprintf ('The result is saved in the files %s \n and %s \n \n', ...
           [output_file '_vel'], [output_file '_press']);
vtk_pts = {linspace(0, 1, 20), linspace(0, 1, 20)};
sp_to_vtk (press, space_p, geometry, vtk_pts, [output_file '_press'], 'press')
sp_to_vtk (vel,   space_v, geometry, vtk_pts, [output_file '_vel'  ], {'velocity', 'divergence'}, {'value', 'divergence'})

% 4.3) PLOT IN MATLAB
[eu, F] = sp_eval (vel, space_v, geometry, vtk_pts);
[X,  Y] = deal (squeeze(F(1,:,:)), squeeze(F(2,:,:)));

figure()
subplot (1,2,2)
eu2 = problem_data.velex (X, Y);
quiver (X, Y, squeeze(eu2(1,:,:)), squeeze(eu2(2,:,:)))
axis equal
title('Exact solution')
subplot (1,2,1)
quiver (X, Y, squeeze(eu(1,:,:)), squeeze(eu(2,:,:)))
axis equal
title('Computed solution')

%!demo
%! ex_stokes_square_sg

%!test
%! problem_data.geo_name = 'geo_square.txt';
%! problem_data.drchlt_sides = 1:4;
%! problem_data.nmnn_sides = [];
%! problem_data.viscosity = @(x, y) ones (size (x));
%! fx = @(x, y) (6 * x + y .* cos(x .* y) + 2 * cos(y) .* sin(x));
%! fy = @(x, y) (x .* cos(x .* y) - 2 * cos(x) .* sin(y));
%! problem_data.f  = @(x, y) cat(1, ...
%!                 reshape (fx (x,y), [1, size(x)]), ...
%!                 reshape (fy (x,y), [1, size(x)]));
%! uxex = @(x, y) (sin(x) .* cos(y));
%! uyex = @(x, y) (-sin(y) .* cos(x));
%! problem_data.h = @(x, y, iside) cat(1, ...
%!                   reshape (uxex (x,y), [1, size(x)]), ...
%!                   reshape (uyex (x,y), [1, size(x)]));
%! problem_data.velex = @(x, y) cat(1, ...
%!                   reshape (uxex (x,y), [1, size(x)]), ...
%!                   reshape (uyex (x,y), [1, size(x)]));
%! problem_data.pressex = @(x, y) (3 * x.^2 + sin(x .* y) - 1.239811742000564725943866);
%! problem_data.gradvelex = @test_stokes_square_bc_graduex;
%! method_data.element_name = 'SG';   % Element type for discretization
%! method_data.degree       = [ 3  3];  % Degree of the splines (pressure space)
%! method_data.regularity   = [ 2  2];  % Regularity of the splines (pressure space)
%! method_data.nsub         = [10 10];  % Number of subdivisions
%! method_data.nquad        = [ 5  5];  % Points for the Gaussian quadrature rule
%! [geometry, msh, space_v, vel, space_p, press] = ...
%!                        solve_stokes (problem_data, method_data);
%! error_l2_p = sp_l2_error (space_p, msh, press, problem_data.pressex);
%! [error_h1_v, error_l2_v] = ...
%!    sp_h1_error (space_v, msh, vel, problem_data.velex, problem_data.gradvelex);
%! assert (msh.nel, 400)
%! assert (space_p.ndof, 169)
%! assert (space_v.ndof, 1152)
%! assert (error_l2_p, 1.95305776265391e-08, 1e-15)
%! assert (error_h1_v, 1.89323099160505e-08, 1e-15)
%! assert (error_l2_v, 2.87241772271343e-10, 1e-15)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [ET_table,EV_table,ETV_index]=edge_tangents_double(V,Ne)
% Edge tangents table
ET_table=zeros(size(V,1)*4,3);
% Edge velocity table
EV_table=zeros(size(V,1)*4,1);
% Edge tangents/velocity index for tables
ETV_index=zeros(size(V,1)*4,2);
ETV_num=0;

% Calculate the tangents and velocity for each edge
Pn=zeros(length(Ne),3); Pnop=zeros(length(Ne),3);
for i=1:size(V,1)
    P=V(i,:);
    Pneig=Ne{i};

    % Find the opposite vertex of each neigbourh vertex.
    % incase of odd number of neigbourhs interpolate the opposite neigbourh
    if(mod(length(Pneig),2)==0)
       for k=1:length(Pneig)
            neg=k+length(Pneig)/2; if(neg>length(Pneig)), neg=neg-length(Pneig); end
            Pn(k,:) = V(Pneig(k),:); Pnop(k,:) = V(Pneig(neg),:);
       end
    else
       for k=1:length(Pneig)
           neg=k+length(Pneig)/2; neg1=floor(neg); neg2=ceil(neg);
           if(neg1>length(Pneig)), neg1=neg1-length(Pneig); end
           if(neg2>length(Pneig)), neg2=neg2-length(Pneig); end
           Pn(k,:) = V(Pneig(k),:); Pnop(k,:) = (V(Pneig(neg1),:)+V(Pneig(neg2),:))/2;
       end
    end

    for j=1:length(Pneig);
        % Calculate length edges of face
        Ec=sqrt(sum((Pn(j,:)-P).^2))+1e-14;
        Eb=sqrt(sum((Pnop(j,:)-P).^2))+1e-14;
        Ea=sqrt(sum((Pn(j,:)-Pnop(j,:)).^2))+1e-14;

        % Calculate face surface area
        s = ((Ea+Eb+Ec)/2);
        h = (2/Ea)*sqrt(s*(s-Ea)*(s-Eb)*(s-Ec))+1e-14;
        x = (Ea^2-Eb^2+Ec^2)/(2*Ea);

        % 2D triangle coordinates
        % corx(1)=0;    cory(1)=0;
        % corx(2)=x;    cory(2)=h;
        % corx(3)=Ea;   cory(3)=0;
        % corx(4)=0;    cory(4)=0;

        % Calculate tangent of 2D triangle
        Np=[-h x]; Np=Np/(sqrt(sum(Np.^2))+1e-14);
        Ns=[h Ea-x]; Ns=Ns/(sqrt(sum(Ns.^2))+1e-14);
        Nb=Np+Ns;
        Tb=[Nb(2) -Nb(1)];

        % Back to 3D coordinates
        Pm=(Pn(j,:)*x+Pnop(j,:)*(Ea-x))/Ea;
        X3=(Pn(j,:)-Pnop(j,:))/Ea;
        Y3=(P-Pm)/h;

        % 2D tangent to 3D tangent
        Tb3D=(X3*Tb(1)+Y3*Tb(2));  Tb3D=Tb3D/(sqrt(sum(Tb3D.^2))+1e-14);

        % Edge Velocity
        Vv=0.5*(Ec+0.5*Ea);

        ETV_num=ETV_num+1;
        ETV_index(ETV_num,:)=[i Pneig(j)];
        ET_table(ETV_num,:)= Tb3D;
        EV_table(ETV_num)=Vv;
    end
end


text:
% Script demonstrating usage of the bpdn function.
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2015-03-05
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'Copyright' and 'License' files
% distributed with the library.

% Signal and dictionary size
N = 512;
M = 4*N;
% Number of non-zero coefficients in generator
L = 32;
% Noise level
sigma = 0.5;

% Construct random dictionary and random sparse coefficients
D = randn(N, M);
x0 = zeros(M, 1);
si = randperm(M);
si = si(1:L);
x0(si) = randn(L, 1);
% Construct reference and noisy signal
s0 = D*x0;
s = s0 + sigma*randn(N,1);

% BPDN for recovery of sparse representation
lambda = 20;
opt = [];
opt.Verbose = 1;
opt.rho = 100;
opt.RelStopTol = 1e-6;
[x1, optinf] = bpdn(D, s, lambda, opt);

figure('position', [100, 100, 1200, 300]);
subplot(1,4,1);
plot(optinf.itstat(:,2));
xlabel('Iterations');
ylabel('Functional value');
subplot(1,4,2);
semilogy(optinf.itstat(:,5));
xlabel('Iterations');
ylabel('Primal residual');
subplot(1,4,3);
semilogy(optinf.itstat(:,6));
xlabel('Iterations');
ylabel('Dual residual');
subplot(1,4,4);
semilogy(optinf.itstat(:,9));
xlabel('Iterations');
ylabel('Penalty parameter');

figure;
plot(x0,'r');
hold on;
plot(x1,'b');
hold off;
legend('Reference', 'Recovered', 'Location', 'SouthEast');
title('Dictionary Coefficients');

% Illustrate restart capability: do 40 iterations and compare with
% 20 iterations and then restart for an additional 20 iterations
opt2 = opt;
opt2.MaxMainIter = 40;
[x2, optinf2] = bpdn(D, s, lambda, opt2);

opt3 = opt2;
opt3.MaxMainIter = 20;
[x3, optinf3] = bpdn(D, s, lambda, opt3);

opt4 = opt3;
opt4.Y0 = optinf3.Y;
opt4.U0 = optinf3.U;
opt4.rho = optinf3.rho;
[x4, optinf4] = bpdn(D, s, lambda, opt4);

figure;
plot(optinf2.itstat(:,1), optinf2.itstat(:,2), 'r');
hold on;
plot(optinf3.itstat(:,1), optinf3.itstat(:,2), 'g');
plot(optinf4.itstat(:,1) + optinf3.itstat(end,1), optinf4.itstat(:,2), 'b');
xlabel('Iterations');
ylabel('Functional value');
legend('Uninterrupted', 'Stop at 20 iterations', 'Restart at 20 iterations');

% Construct ensemble of L-sparse coefficients and corresponding signals
K = 8;
X0 = zeros(M, K);
for l = 1:K,
  si = randperm(M);
  si = si(1:L);
  X0(si,l) = randn(L, 1);
end
S0 = D*X0;
S = S0 + sigma*randn(N, K);

% BPDN for simultaneous recovery of sparse coefficient matrix
lambda = 20;
opt = [];
opt.Verbose = 1;
opt.rho = 100;
opt.RelStopTol = 1e-6;
[X1, optinf] = bpdn(D, S, lambda, opt);

figure;
subplot(1,2,1);
imagesc(X0);
title('Reference');
subplot(1,2,2);
imagesc(X1);
title('Recovered');

author: thfylsty


text:
%% Ellipse Based Shape Parameters
%
%%
% In this section we discuss geometric properties of grains that are
% related to ellipses fitted to the grains. Additionally to the orientation
% |omega|, and the lengths |a|, |b| of the long axis and short axes that
% are computed by the command <grain2d.fitEllipse.html |[omega,a,b] =
% grains.fitEllipse|>  the following properties based on the fitted
% ellipses are avaiable.
%
% || <grain2d.longAxis.html |longAxis|> || long axis as @vector3d || <grain2d.shortAxis.html |shortAxis|>  || short axis as @vector3d ||
% || <grain2d.centroid.html |centroid|> || midpoint  || <grain2d.aspectRatio.html |aspectRatio|>  || long axis / short axis ||
%
% In order to demonstrate these properties we start by reconstructing the
% grain structure from a sample EBSD data set.

% load sample EBSD data set
mtexdata forsterite silent

% reconstruct grains and smooth them
[grains, ebsd.grainId] = calcGrains(ebsd('indexed'),'angle',5*degree);
ebsd(grains(grains.grainSize<10)) = [];
[grains, ebsd.grainId] = calcGrains(ebsd('indexed'),'angle',5*degree);
grains(grains.isBoundary) = [];

grains=smooth(grains('indexed'),10,'moveTriplePoints');

% plot the grains
plot(grains,'micronbar','off','lineWidth',2)

%% Fit Ellipses
%
% The basic command for fitting ellipses is <grain2d.fitEllipse
% |fitEllipse|>

[omega,a,b] = grains.fitEllipse;

plotEllipse(grains.centroid,a,b,omega,'lineColor','w','linewidth',2)

%%
% The returned variable |omega| is the angle describing the rotation of the
% ellipses and |a| and |b| are the length of the longest and shortest half
% axis. The midpoints of the ellipses can be computed by the command
% <grain2d.centroid.html |grains.centroid|>.  Note, that the ellipses are
% scaled such that the area of the ellipse coincides with the actual grain
% area. Alternatively, one can also scale the ellipse to fit the boundary
% length by using the option |boundary|.
%
%% Long and Short Axes
%
% The direction of the long and the short axis of the fitted ellipse can be
% obtained by the comands <grain2d.longAxis.html |grains.longAxis|> and
% <grain2d.shortAxis.html |grains.shortAxis|>. These directions are only
% well defined if the fitted ellipse is not to close to a perfect circle. A
% measure for how distinct the ellipse is from a perfect circle is the
% <grain2d.aspectRatio.html aspect ratio> which is defined as the quotient
% $a/b$ between the longest and the shortest axis. For a perfect circle
% the apect ratio is $1$ and increases to infinity when the ellipse becomes
% more and more elongated.
%
% Lets colorize the grains by their apect ratio and plot on top the long
% axis directions:

% visualize the aspect ratio
plot(grains,grains.aspectRatio,'linewidth',2,'micronbar','off')
setColorRange([0,4])
mtexColorbar('title','aspect ratio')

% and on top the long axes
hold on
quiver(grains,grains.longAxis,'Color','white')
hold off

%% Shape perfered orientation
%
% If we look at grains, we might wonder if there is a characteristic
% difference in the grain shape fabric between e.g. Forsterite and
% Enstatite. In contrast to crystal prefered orientations which which
% describe on the alignment of the atome lattices the shape prefered
% orientation (SPO) describes the algnment of the grains by shape in the
% bulk fabric.
%
% *Long Axis Distribution*
%
% The most direct way to analyse shape prefered orientations are rose
% diagrams of the distribution of the grain long axes. For those diagrams
% it is useful to weight the long axis by the grain area such that larger
% grains have a bigger impact on the distribution and by the aspect ratio
% as for grains with a small aspect ratio the long axis is not so well
% defined.

numBin = 50;

subplot(1,2,1)
weights = grains('forsterite').area .* (grains('forsterite').aspectRatio-1);
histogram(grains('forsterite').longAxis,numBin, 'weights', weights)
title('Forsterite')

subplot(1,2,2)
weights = grains('enstatite').area .* (grains('enstatite').aspectRatio - 1);
histogram(grains('enstatite').longAxis,numBin,'weights',weights)
title('Enstatite')

%%
% Instead of the histogram we may also fit a circular density distribution
% to the to the long axes using the command <calcDensity.thml
% |calcDensity|>.

tdfForsterite = calcDensity(grains('forsterite').longAxis,...
  'weights',norm(grains('forsterite').longAxis),'halfwidth');

tdfEnstatite = calcDensity(grains('enstatite').longAxis,...
  'weights',norm(grains('enstatite').longAxis));

plotSection(tdfForsterite, vector3d.Z, 'linewidth', 3)

hold on
plotSection(tdfEnstatite, vector3d.Z, 'linewidth', 3)
hold off

%%
%

close all
[freq,bc] = calcTDF(grains('fo'),'binwidth',3*degree);
plotTDF(bc,freq/sum(freq));

[freq,bc] = calcTDF(grains('en'),'binwidth',3*degree);
hold on
plotTDF(bc,freq/sum(freq));
hold off
legend('Forsterite','Enstatite','Location','eastoutside')
mtexTitle('long axes')

%% *Shortest Caliper Distribution*
%
% Alternatively, we may wonder if the common long axis of grains is does
% suitably represented by the direction normal to the shortest caliper of
% the grains. This can particularly be the case for aligned rectangular
% particles. The command <calcTDF.html |calcTDF|> also takes a list of
% angles and a list of weights or lengths as input

cPerpF = caliper(grains('fo'),'shortestPerp');
cPerpE = caliper(grains('en'),'shortestPerp');

[freqF,bcF] = calcTDF(cPerpF.rho, 'weights',cPerpF.norm, 'binwidth',3*degree);
[freqE,bcE] = calcTDF(cPerpE.rho, 'weights',cPerpE.norm, 'binwidth',3*degree);

plotTDF(bcF,freqF/sum(freqF));
hold on
plotTDF(bcE,freqE/sum(freqE));
hold off
legend('Forsterite','Enstatite','Location','eastoutside')

%%
% We can also smooth the functions with a wrapped Gaussian

pdfF = circdensity(bcF, freqF, 5*degree,'sum');
pdfE = circdensity(bcE, freqE, 5*degree,'sum');

plotTDF(bcF,pdfF);
hold on
plotTDF(bcE,pdfE);
hold off
mtexTitle('n.t.s. density estimate')
legend('Forsterite','Enstatite','Location','eastoutside')

%%
% Because best fit ellipses are always symmetric and the projection
% function of an entire grain always only consider the convex hull, grain
% shape fabrics can also be characterized by the the length weighted rose
% diagram of the directions of grain boundary segments.

[freqF,bcF] = calcTDF(grains('fo').boundary);
plotTDF(bcF,freqF/sum(freqF));
pdfF = circdensity(bcF, freqF, 5*degree,'sum');
hold on
plotTDF(bcF,pdfF);
hold off
mtexTitle('Forsterite grain boundaries')
nextAxis
[freqE,bcE] = calcTDF(grains('en').boundary);
plotTDF(bcE,freqE/sum(freqE));
pdfE = circdensity(bcE, freqE, 5*degree,'sum');
hold on
plotTDF(bcE,pdfE);
hold off
mtexTitle('Enstatite grain boundaries')

%% Characteristic Shape
%
% Note that this distribution is very prone to inherit artifacts based on
% the fact that most EBSD maps are sampled on a regular grid. We tried to
% overcome this problem by heavily smoothing the grain boundary. The little
% peaks at 0 and 90 degree are very likely still related to this sampling
% artifact.
%
% If we just add up all the individual elements of the rose diagram in
% order of increasing angles, we derive the characteristic shape. It can be
% regarded as to represent the average grain shape.

[csAngleF, csRadiusF] = characteristicShape(bcF,freqF);
[csAngleE, csRadiusE] = characteristicShape(bcE,freqE);

close all
plotTDF(csAngleF,csRadiusF,'nolabels');
hold on
plotTDF(csAngleE,csRadiusE,'nolabels');
hold off
legend('Forsterite','Enstatite','Location','eastoutside')

%%
% We may wonder if these results are significantly different or not
% TODO: get deviation from an ellipse etc

author: mtex-toolbox

repo: mtex

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


text:
function state = dss_core_mim(state)
% MIM DSS algorithm
%   state = dss_core_mim(state)
%     Calculates linear denoising source separation with MIM.
%     Performs calculation defined by state structure and
%     returns the result as state structure.
%     Used CPU time is recorded in the state.

% Copyright (C) 2004, 2005 DSS MATLAB package team (dss@cis.hut.fi).
% Distributed by Laboratory of Computer and Information Science,
% Helsinki University of Technology. http://www.cis.hut.fi/projects/dss/.
% $Id: dss_core_pca.m,v 1.13 2005/12/02 12:23:18 jaakkos Exp $

dss_message(state,2,'Extracting components in PCA DSS\n');

% -- Initialize local variables
start_time = cputime;
interrupt_iteration = 0;
interrupt_component = 0;
% Dimension of the whitened data
if iscell(state.Y)
  wdim = size(state.Y{1},1);
else
  wdim = size(state.Y, 1);
end
sdim = state.sdim;

% -- Denoising
[state.denf.params, Y] = feval(state.denf.h, state.denf.params, state.Y, state);

% -- PCA
if iscell(Y)
  C = cellcov(Y,2,0);
  %[E,D] = eig(cellcov(Y,2,0));
else
  C = (Y*Y')/size(Y,2);
  %[E,D] = eig(Y * Y' / size(Y,2));
end

i1 = state.indx==1;
i2 = state.indx==2;

% create whitening matrix T
T = C;
T(i1, i2) = 0;
T(i2, i1) = 0;
T = inv(sqrtm(real(T)));

E = imag(T(i1, i1) * C(i1, i2) *T(i2, i2)');

[e1, d1] = eig(E*E');
[e2, d2] = eig(E'*E);

d1 = diag(d1);
d2 = diag(d2);
n  = min(sum(i1), sum(i2));

% sort eigenvalues in descending order
[d1,order] = sort(-d1);
d1 = -d1(1:n);
e1 = e1(:,order(1:n));

[d2,order] = sort(-d2);
d2 = -d2(1:n);
e2 = e2(:,order(1:n));

state.D = [d1;d2];
state.sdim = numel(state.D);
state.W = blkdiag(e1,e2)'/T;
state.S = state.W * state.Y;

% -- record the used cpu time
if ~isfield(state, 'cputime'); state.cputime = 0; end
state.cputime = state.cputime + cputime - start_time;


text:
function [X_norm, mu, sigma] = featureNormalize(X)
%FEATURENORMALIZE Normalizes the features in X
%   FEATURENORMALIZE(X) returns a normalized version of X where
%   the mean value of each feature is 0 and the standard deviation
%   is 1. This is often a good preprocessing step to do when
%   working with learning algorithms.

mu = mean(X);
X_norm = bsxfun(@minus, X, mu);

sigma = std(X_norm);
X_norm = bsxfun(@rdivide, X_norm, sigma);

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

end

author: Ayatans

repo: Machine-Learning-homework


text:
%% EXAMPLE 2: Plot SPOD spectrum and inspect SPOD modes.
%  The large-eddy simulation data provided along with this example is a
%  subset of the database of a Mach 0.9 turbulent jet described in [1] and
%  was calculated using the unstructured flow solver Charles developed at
%  Cascade Technologies. If you are using the database in your research or
%  teaching, please include explicit mention of Brès et al. [1]. The test
%  database consists of 5000 snapshots of the symmetric component (m=0) of
%  a round turbulent jet. A physical interpretaion of the SPOD results is
%  given in [2], and a comprehensive discussion and derivation of SPOD and
%  many of its properties can be found in [3].
%
%   References:
%     [1] G. A. Brès, P. Jordan, M. Le Rallic, V. Jaunet, A. V. G.
%         Cavalieri, A. Towne, S. K. Lele, T. Colonius, O. T. Schmidt,
%         Importance of the nozzle-exit boundary-layer state in subsonic
%         turbulent jets, J. of Fluid Mech. 851, 83-124, 2018
%     [2] Schmidt, O. T. and Towne, A. and Rigas, G. and Colonius, T. and
%         Bres, G. A., Spectral analysis of jet turbulence, J. of Fluid Mech. 855, 953–982, 2018
%     [3] Towne, A. and Schmidt, O. T. and Colonius, T., Spectral proper
%         orthogonal decomposition and its relationship to dynamic mode
%         decomposition and resolvent analysis, J. of Fluid Mech. 847, 821–867, 2018
%
% O. T. Schmidt (oschmidt@ucsd.edu), A. Towne, T. Colonius
% Last revision: 5-Sep-2022 (OTS)

clc, clear variables
addpath('utils')
load(fullfile('jet_data','jetLES.mat'),'p','x','r','dt');

%% SPOD of the test database.
%   Calculate the SPOD of the data matrix 'p' and use the timestep 'dt'
%   between snapshots to obtain the physical frequency 'f'. 'L' is the
%   matrix of modal energies, as before, and 'P' the data matrix of SPOD
%   modes. We leave all other options empty for now.
[L,P,f] = spod(p,[],[],[],dt);

%   First, we plot the SPOD spectrum.
figure
loglog(f,L)
xlabel('frequency'), ylabel('SPOD mode energy')

%   Second, we visualize the 1st and 2nd SPOD modes at three frequencies.
figure
count = 1;
for fi = [10 15 20]
    for mi = [1 2]
        subplot(3,2,count)
        contourf(x,r,real(squeeze(P(fi,:,:,mi))),11,'edgecolor','none'), axis equal tight, caxis(max(abs(caxis))*[-1 1])
        xlabel('x'), ylabel('r'), title(['f=' num2str(f(fi),'%.2f') ', mode ' num2str(mi) ', \lambda=' num2str(L(fi,mi),'%.2g')])
        xlim([0 10]); ylim([0 2])
        count = count + 1;
    end
end

%% Animate the same modes.
%   Note how all wavepackets travel at approximately the same phase
%   speed c_ph. The reason is that their streamwise wavenumber k_x changes
%   with frequency such that c_ph = omega/k_x is approximately constant.
figure
nt      = 30;
T       = 1/f(10);              % period of the 10th frequency
time    = linspace(0,T,nt);     % animate over one period
count = 1;
for ti = 1:nt
    for fi = [10 15 20]
        for mi = [1 2]
            subplot(3,2,count)
            pcolor(x,r,real(squeeze(P(fi,:,:,mi)*exp(2i*pi*f(fi)*time(ti))))), shading interp, axis equal tight, caxis(max(abs(caxis))*[-1 1])
            xlabel('x'), ylabel('r'), title(['f=' num2str(f(fi),'%.2f') ', mode ' num2str(mi) ', \lambda=' num2str(L(fi,mi),'%.2g')])
            xlim([0 10]); ylim([0 2])
            count = count + 1;
            hold on
        end
    end
    drawnow
    hold off
    count = 1;
end

author: SpectralPOD


text:
function grid_weight = sparse_grid_weights_ofn ( dim_num, level_max, rule, ...
  point_num, grid_index )

%*****************************************************************************80
%
%% SPARSE_GRID_WEIGHTS_OFN computes sparse grid weights based on a OFN 1D rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Input, integer RULE, the index of the rule.
%    1, "CC", Clenshaw Curtis Closed Fully Nested rule.
%    2, "F1", Fejer 1 Open Fully Nested rule.
%    3, "F2", Fejer 2 Open Fully Nested rule.
%    4, "GP", Gauss Patterson Open Fully Nested rule.
%    5, "GL", Gauss Legendre Open Weakly Nested rule.
%    6, "GH", Gauss Hermite Open Weakly Nested rule.
%    7, "LG", Gauss Laguerre Open Non Nested rule.
%
%    Input, integer POINT_NUM, the total number of points in
%    the grids.
%
%    Input, integer GRID_INDEX(DIM_NUM,POINT_NUM), a list of
%    point indices, representing a subset of the product grid of level
%    LEVEL_MAX, representing (exactly once) each point that will show up in a
%    sparse grid of level LEVEL_MAX.
%
%    Output, real GRID_WEIGHT(POINT_NUM), the weights
%    associated with the sparse grid points.
%
  if ( level_max == 0 )
    grid_weight(1:point_num) = 2.0^dim_num;
    return
  end

  grid_weight(1:point_num) = 0.0;

  level_min = max ( 0, level_max + 1 - dim_num );

  for level = level_min : level_max
%
%  The middle loop generates the next partition LEVEL_1D(1:DIM_NUM)
%  that adds up to LEVEL.
%
    level_1d = [];
    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );
%
%  Transform each 1D level to a corresponding 1D order.
%
      order_1d = level_to_order_open ( dim_num, level_1d );
%
%  The product of the 1D orders gives us the number of points in this grid.
%
      order_nd = prod ( order_1d(1:dim_num) );
%
%  Generate the indices of the points corresponding to the grid.
%
      grid_index2 = multigrid_index_ofn ( dim_num, order_1d, order_nd );
%
%  Compute the weights for this grid.
%
      grid_weight2 = product_weights ( dim_num, order_1d, order_nd, rule );
%
%  Adjust the grid indices to reflect LEVEL_MAX.
%
      grid_index2 = multigrid_scale_open ( dim_num, order_nd, level_max, ...
        level_1d, grid_index2 );
%
%  Now determine the coefficient.
%
      coeff = r8_mop ( level_max - level ) ...
        * r8_choose ( dim_num - 1, level_max - level );

      for point2 = 1 : order_nd

        for point = 1 : point_num

          if ( all ( ...
            grid_index2(1:dim_num,point2) == grid_index(1:dim_num,point) ...
          ) )
            grid_weight(point) = grid_weight(point) ...
              + coeff * grid_weight2(point2);
            break
          end

        end

      end

      if ( ~more )
        break
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function H = lbfgsinit(n,k,dscale)
%lbfgsinit  initializes an L-BFGS matrix.
%
%   H = lbfgsinit(n,k,dscale)
%
%   See also lbfgsadd, lbfgsdel, lbfgshprod, lbfgsupdate.

% $Id: lbfgsinit.m 385 2008-09-02 22:22:41Z mpf $

if nargin < 3 || isempty(dscale)
   dscale = 1;
end

H.jNew  = 1;
H.jOld  = 1;
H.jMax  = k;
H.S     = zeros(n,k);
H.Y     = zeros(n,k);

% Data specific to H
H.gamma = dscale;
H.r     = zeros(1,k);

% Data specific to B
H.delta = 1 / dscale;
H.valid = false(k,1);
H.rank  = 0;
H.STS   = zeros(k,k);
H.L     = zeros(k,k);
H.D     = zeros(k,k);
H.M     = zeros(2*k,2*k);
H.ML    = [];
H.MU    = [];

% Status information and options
H.status = 0; % 0=updated, 1=corrected, 2=no update

author: mpf


text:
function [kN] = kip2kN(kip)
% Convert force from kip to kilonewtons.
% Chad A. Greene 2012
kN = kip*4.4482216;

author: Sable

repo: mcbench-benchmarks


text:
function M = convectionUpwindTermCylindrical1D(u)
% This function uses the upwind scheme to discretize a 1D
% convection term in the form \grad (u \phi) where u is a face vactor
% It also returns the x and y parts of the matrix of coefficient.
%
% SYNOPSIS:
%   M = convectionUpwindTermCylindrical1D(u)
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% extract data from the mesh structure
Nr = u.domain.dims(1);
G = 1:Nr+2;
DXp = u.domain.cellsize.x(2:end-1);
rp = u.domain.cellcenters.x;
rf = u.domain.facecenters.x;

% define the vectors to store the sparse matrix data
iix = zeros(3*(Nr+2),1);
jjx = zeros(3*(Nr+2),1);
sx = zeros(3*(Nr+2),1);

% extract the velocity data
% note: size(ux) = [1:m+1, 1:n] and size(uy) = [1:m, 1:n+1]
ux = u.xvalue;

% reassign the east, west, north, and south velocity vectors for the
% code readability
ue = ux(2:Nr+1);		uw = ux(1:Nr);
re = rf(2:Nr+1);     rw = rf(1:Nr);

% find the velocity direction for the upwind scheme
ue_min = min(ue,0);	ue_max = max(ue,0);
uw_min = min(uw,0);	uw_max = max(uw,0);

% calculate the coefficients for the internal cells
AE = reshape(re.*ue_min./(DXp.*rp),Nr,1);
AW = reshape(-rw.*uw_max./(DXp.*rp),Nr,1);
APx = reshape((re.*ue_max-rw.*uw_min)./(DXp.*rp),Nr,1);

% correct for the cells next to the boundary
% Left boundary:
APx(1) = APx(1)-rw(1)*uw_max(1)/(2*rp(1)*DXp(1));   AW(1) = AW(1)/2;
% Right boundary:
AE(end) = AE(end)/2;    APx(end) = APx(end)+re(end)*ue_min(end)/(2*rp(end)*DXp(end));

% build the sparse matrix based on the numbering system
rowx_index = reshape(G(2:Nr+1),Nr,1); % main diagonal x
iix(1:3*Nr) = repmat(rowx_index,3,1);
jjx(1:3*Nr) = [reshape(G(1:Nr),Nr,1); reshape(G(2:Nr+1),Nr,1); reshape(G(3:Nr+2),Nr,1)];
sx(1:3*Nr) = [AW; APx; AE];

% build the sparse matrix
kx = 3*Nr;
M = sparse(iix(1:kx), jjx(1:kx), sx(1:kx), Nr+2, Nr+2);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function img = mrAnatMontage(img, xform, acpcSlices, fname, figHandle)
% montageRgbImg = mrAnatMontage(img, xform, acpcSlices, fname, figHandle)
%
% Quick hack for displaying nice anatomy montages with ac-pc slice
% labels. Must be a volume, but can be an intensity image (will be
% displayed as grayscale) or an [x,y,z,3] rgb volume.
%
% HISTORY:
% 2006.08.08 RFD: wrote it.
%

if(~exist('fname','var')) fname = ''; end
if(~exist('figHandle','var')) figHandle = figure; end

if(size(img,4)==1)
  img = repmat(img,[1,1,1,3]);
end
% reorient so that the eyes point up
img = flipdim(permute(img,[2 1 3 4]),1);
for(ii=1:length(acpcSlices)) slLabel{ii} = sprintf('Z = %d',acpcSlices(ii)); end
[t,r,s] = affineDecompose(xform);
slImg = inv(xform)*[zeros(length(acpcSlices),2) acpcSlices' ones(length(acpcSlices),1)]';
slImg = round(slImg(3,:));
img = makeMontage3(img, slImg, s(1), 0, slLabel,[],figHandle);
if(~isempty(fname)) mrUtilPrintFigure(fname); end

% be kind to those who forget a semicolon:
if(nargout<1) clear img; end
return;


text:
function order = IDX2order( IDX )
% Converts class labels into an ordering.
%
% Creates an ordering order such that IDX(order)=[1 1...1 2...2 ... k...k].
% All points within a class retain the ordering in which they originally
% appeared.  Also, Xb = X(order,:) has cluster labels IDX(order), ie
% adjacent elements in X typically belong to the same cluster.
%
% USAGE
%  order = IDX2order( IDX )
%
% INPUTS
%  IDX     - cluster membership [see kmeans2.m]
%
% OUTPUTS
%  order   - n-by-1 vector containing a new ordering for the points.
%
% EXAMPLE
%  order = IDX2order( [1 1 3 1 2 2] )  % should be: [1 2 4 5 6 3]
%
% See also DISTMATRIXSHOW
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

k = max(IDX);  n = length(IDX);
order = zeros(1,n);  count = 0;
for i=1:k
  locs = (IDX==i); orderi = cumsum(locs);
  order(locs) = orderi(locs) + count;
  count = count+sum(locs);
end
[dis,order] = sort(order);

author: zouchuhang

repo: LayoutNet


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  June 2013
% ------------------------------------------------------------------------

%% Dummy dummy
lp = [ 1  1  1
       2  3  2
       2  2  2];
ms_matrix = [ 1  2  4
              3  4  5];

 cands     = [1 2];
[cands_hf] = hole_filling(lp, ms_matrix, cands);
assert(isequal(cands_hf,[1 2 3]))

%% Dummy
lp = [ 1  2  3  4  5
       6  7  8  9 10
      11 12 13 14 15
      16 17 18 19 20];

    % 21 21 21 21 24
    % 22 25 25 21 24
    % 22 26 26 23 24
    % 23 23 23 23 24
ms_matrix = [ 1  2  3  4  9 21
              6 11  0  0  0 22
             16 17 18 19 14 23
             10 15  5 20  0 24
              7  8  0  0  0 25
             12 13  0  0  0 26
             22 23 24  0  0 27
             26 27  0  0  0 28
             21 25  0  0  0 29
             28 29  0  0  0 30];

 cands     = [21 22 23];
[cands_hf, cands_comp] = hole_filling(lp, ms_matrix, cands);
assert(isequal(cands_hf  ,[21 22 23 25 26]))
assert(isequal(cands_comp,24))

 cands2     = [21 22 23 24 25];
[cands_hf, cands_comp] = hole_filling(lp, ms_matrix, cands2);
assert(isequal(cands_hf  ,[21 22 23 24 25 26]))
assert(isequal(cands_comp,30))

 cands3    = [1 2 3];
[cands_hf, cands_comp] = hole_filling(lp, ms_matrix, cands3);
assert(isequal(cands_hf  ,[1 2 3]))
assert(isequal(cands_comp,[4 9 25 28]))

%% Real test
im_id = '2008_000009';
load(fullfile(root_dir,'datasets','pascal2012','gPb_mUCM','multi', [im_id '.mat']));

curr_hier = ucm2hier(ucm2);
ths{1}.start_ths = curr_hier.start_ths';
ths{1}.end_ths = curr_hier.end_ths';
ms{1} = curr_hier.ms_matrix;
lps = curr_hier.leaves_part;

[f_lp,f_ms,cands] = full_cands_from_hiers(lps,ms,ths,[500 500 500 500]');

tic
[cands_hf, cands_comp] = hole_filling(double(f_lp), double(f_ms), cands);
toc

tic
% Get masks
masks    = cands2masks(cands, f_lp, f_ms);
masks_hf2= false(size(masks));

% Perform hole filling by morphology
for ii=1:size(masks,3)
    masks_hf2(:,:,ii) = (imfill(masks(:,:,ii),'holes')>0);
end
toc

% Is equal?
masks_hf = cands2masks(cands_hf, f_lp, f_ms);
assert(isequal(masks_hf,masks_hf2))

% Was there any hole?
assert(~isequal(masks_hf,masks))

% Check complementaries
n_regs = f_ms(end,end);
masks_comp = cands2masks(cands_comp, f_lp, f_ms);
for ii=1:size(masks,3)
    if isequal(cands_comp(ii),n_regs)
        assert(unique(masks_hf(:,:,ii))==1)
    else
        tmp = double(masks_comp(:,:,ii))+double(masks_hf(:,:,ii));
        assert(unique(tmp)==1)
    end
end

%% Show some results
% id = 2000;
% figure;
% subplot(1,3,1);imshow(masks(:,:,id)>0)
% subplot(1,3,2);imshow(masks_hf(:,:,id)>0)
% subplot(1,3,3);imshow(masks_comp(:,:,id)>0)

author: s-gupta

repo: rcnn-depth

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


text:
report_this_filefun(mfilename('fullpath'));

%l = isnan(tmap);
%tmap(l) = 1;

l = tmap < 100;
tmap(l) = 100;

[lat,lon] = meshgrat(tmap,tmapleg);
%[smap,smapleg] = country2mtx('switzerland',100);
%[lat0, lon0] = meshgrat(smap,smapleg);

% tmap = km2deg(tmap/1);
[X , Y]  = meshgrid(gx,gy);

%sw = interp2(lon0,lat0,smap,lon,lat);

re4 = re3;

ren = interp2(X,Y,re4,lon,lat);

mi = min(min(ren));
l =  isnan(ren);
ren(l) = mi-20;

figure_w_normalized_uicontrolunits('pos',[150 500 1000 700])

hold on; axis off
axesm('MapProjection','eqaconic','MapParallels',[],...
    'MapLatLimit',[33.75 35.03],'MapLonLimit',[-116.8 -115.9])

meshm(ren,tmapleg,size(tmap),tmap);

daspectm('m',4);
tightmap
%view([10 30])
hl = camlight ; lighting phong
set(gca,'projection','perspective');

aa = a;

j = jet;
%j = j(64:-1:1,:);
j = [ [ 0.85 0.9 0.9] ; j];
caxis([ min(min(re4)) max(max(re4)) ]);

colormap(j);
caxis([0.60 1.5]);
axis off; set(gcf,'color','w')
ax = axis;
axis([ax(1:4) 0 5000])

setm(gca,'ffacecolor','w')
setm(gca,'fedgecolor','k','flinewidth',3);

setm(gca,'mlabellocation',0.5)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',0.5)
setm(gca,'parallellabel','on','Labelrotation','on')
setm(gca,'Fontcolor','k','Fontweight','bold','FontSize',14,'Labelunits','dm')

h5 = colorbar;
set(h5,'position',[0.65 .13 0.015 0.3],'TickDir','out','Ycolor','k','Xcolor','k',...
    'Fontweight','bold','FontSize',14,'Ticklength',[0.025 0.02]);
set(gcf,'Inverthardcopy','off');

a = aa;

author: CelsoReyes


text:
function generator_order1_norad_no(filenameout,GeneratedYears,TempScheme,idistr,...
    ap00,ap10,alambda,nu,A,B,aC0,aC1,aC2,aD1,aD2,sC0,sC1,sC2,sD1,...
    sD2,PrecipThreshold,MarkovChainOrder)
%
%  THIS PROGRAM GENERATES DAILY TIME SERIES OF PRECIPITATION, MAX AND MIN
%  TEMPERATURES. PARAMETERS OF THE MODEL HAVE BEEN CALCULATED BY THE
%  WEATHER ANALYZER PROGRAM. A RANDOM GENERATOR APPROACH IS USED TO CREATE
%  SYNTHETIC TIME SERIES

% the preciptiation occurrence is generated using the first order markov
% chain and the parameters are not smoothed from 14 days to daily scale with
% this weather generator

% 'filenameinter' is a string and is the filename that contains all variables
% saved from the analyzer (analyzer_order1_norad.m).
%
% this version uses a first order Markov chain for precipitation occurence
%
% 'filenameout' is the filename (string) to save the generated variables
% gP, gTmax, gTmin

Y=GeneratedYears;
h = waitbar(0,'Initialisation...');

A=A';
B=B';
waitbar(0.1,h);
%
% generation of time series

% time series for precipitation
%
%  create vectors of daily p00, p10, lbd values from Fourier series
%
day=1:365;
tot=365*Y;
t=1:tot;
T=365/(2*pi);
waitbar(0.2,h);
p00=zeros(Y,356);
p10=zeros(Y,365);
mm=1;
nn=1;
for i=1:14:364
    p00(1,mm:mm+13)=ap00(nn);
    p10(1,mm:mm+13)=ap10(nn);
    mm=mm+14;
    nn=nn+1;
end
p00(1,365)=ap00(26);
p10(1,365)=ap10(26);

for i=2:Y
    p00(i,:)=p00(1,:);
    p10(i,:)=p10(1,:);
end
p00=p00';
p00=reshape(p00,1,[]);

p10=p10';
p10=reshape(p10,1,[]);

waitbar(0.5,h);
p01=1-p00;

waitbar(0.8,h);
p11=1-p10;
%
% add one element to vectors p00, p01, p10 and p11 to avoid putting an if
% statement in the loop to generate time series of wet dry days
%
p00(tot+1)=0;
p01(tot+1)=p00(tot+1);
p10(tot+1)=p00(tot+1);
p11(tot+1)=p00(tot+1);

lbd=zeros(Y,356);
mm=1;
nn=1;
for i=1:14:364
    lbd(1,mm:mm+13)=alambda(nn);
    mm=mm+14;
    nn=nn+1;
end
lbd(1,365)=alambda(26);

for i=2:Y
    lbd(i,:)=lbd(1,:);
end
lbd=lbd';
lbd=reshape(lbd,1,[]);

if idistr == 2
    anu=nu;
    nu=zeros(Y,356);
    mm=1;
    nn=1;
    for i=1:14:364
        nu(1,mm:mm+13)=anu(nn);
        mm=mm+14;
        nn=nn+1;
    end
    nu(1,365)=anu(26);
    for i=2:Y
        nu(i,:)=nu(1,:);
    end
    nu=nu';
    nu=reshape(nu,1,[]);
end
waitbar(1,h);
close(h);
%
% generate time series of dry (X=0) and wet (X=1) days
% pn is the marginal distribution of X on day n
%
% assume that day 0 is dry.  state=0 when day is dry,  state=1 when day is wet
%
state=0;
%
% loop to generate time series
%
n=Y*365;
% iwet=0; % fbri7
X=zeros(1,n); % fbri7

h = waitbar(0,'Generation of precipitation occurrence...');
wbt=round(n/100);

for i=1:n
   % generate random number from uniform distribution
    Ru(i,1)=rand(1);
	% establish which of p00 or p10 tu use
    if state == 0
        % day i-1 is dry, use p00
        Pr=p00(i);
    else
        % day i-1 is wet, use p10
        Pr=p10(i);
    end
    % establish if day i is dry or wet
   if Ru(i,1) <= Pr;
      % dry day
      state=0;
   else
      % wet day
      X(i)=1;
      state=1;
   end
   if i/wbt - round(i/wbt) == 0   % update waitbar at every 100th interval
       waitbar(i/n,h);
   end
end

close(h);
clear p00 p01 p10 p11;

XX=find(X==1);
nw=length(XX);
%
% generate time series of precipitation amounts
%
tol=0.00001;    % tolerance criteria for iteration
r=zeros(1,n);

h = waitbar(0,'Generation of precipitation quantity...');
wbt=round(nw/100);

if idistr == 1
    %
    % assuming exponential distribution
    %
    for i=1:nw
       P=rand(1);
       r(XX(i))=-log(1-P)/lbd(XX(i));
       if i/wbt - round(i/wbt) == 0   % update waitbar every 1/100th
            waitbar(i/(nw*1.1),h);  % waitbar at 90% after this step
       end
    end
else
    %
    % assuming 2-parameter gamma distribution
    %
    for i=1:nw  % for each wet day, draw a random number and get precip using
                % inverse of gamma CDF
       % Px=rand(1);
       Px(i,1)=rand(1);
       r(XX(i))=gaminv(Px(i,1),nu(XX(i)),1/lbd(XX(i)));
       if i/wbt - round(i/wbt) == 0   % update waitbar every 1/100th
           waitbar(i/(nw*1.1),h);      % waitbar at 90% after this step
       end
    end
end

clear lbd;
%
%  time series of min temperature and max temperature
%
%  we first generate the fourier estimates of average and standard deviations
%  of the time series
%
for j=1:4
   ay(j,:)=aC0(j)+aC1(j)*sin(t/T+aD1(j))+aC2(j)*sin(2*t/T+aD2(j));
   sy(j,:)=sC0(j)+sC1(j)*sin(t/T+sD1(j))+sC2(j)*sin(2*t/T+sD2(j));
   waitbar(0.9+j*0.025,h);
end
clear t;
%
%  procedure is started by assuming that the residuals are equal to 0
%  random component eps is normally distributed N(0,1)
%
res=[0; 0];
ksi=zeros(2,n);  % fbri7

close(h);
h = waitbar(0,'Generation of Tmin, Tmax...');
wbt=round(n/100);

for i=1:n
   eps=randn(2,1);
   res=A*res+B*eps;
   ksi(:,i)=res;
   if i/wbt - round(i/wbt) == 0   % update waitbar every 1/100th
       waitbar(i/(n*1.33),h);          % waitbar at 75% when this step is over
   end
end
%
% the means and standard deviations obtained by Fourier time series are conditioned
% on the wet or dry status of the day determined by using the Markov chain model
%
v=ones(2,1);
XX=kron(v,X);
waitbar(0.78,h);
clear X;

cay=XX.*ay(1:2,:)+(1-XX).*ay(3:4,:);
waitbar(0.83,h);
clear ay;
csy=XX.*sy(1:2,:)+(1-XX).*sy(3:4,:);
waitbar(0.86,h);
clear sy XX;
%
if TempScheme==1;
    % the daily values of the three weather variables are found by multiplying the
    % residuals by the standard deviation and adding the mean
    %
    Xp=ksi.*csy+cay;
    clear cay csy ksi;
    Tmax=Xp(1,:);
    Tmin=Xp(2,:);
    clear Xp;
    % insure that Tmin is always smaller than Tmax (diff arbitrarily set up at 1)
    t0=find(Tmax-Tmin < 1);
    Tmin(t0)=Tmax(t0)-1;
elseif TempScheme==2;
    % the Tmax and Tmin are generated conditioned on each other (Jie Chen modified)
    % the smaller standard deviation of Tmax or Tmin is used as a base, and the
    % other parameter is generated conditioned on the chosen parameter. If the
    % standard deviation of Tmax is larger than or equal to the standard
    % deviation of Tmin, daily temperatures are generated by:
    % Tmin=Mean(min)+Std(min)*rand
    % Tmax=Tmin+(Mean(max)-Mean(min))+(Std(max)^2-Std(min)^2)^0.5*rand
    % If the standard deviation of Tmax is less than those of Tmin, daily
    % temperatures are genareted by:
    % Tmax=Mean(max)+Std(max)*rand
    % Tmin=Tmax-(Mean(max)-Mean(min))-(Std(min)^2-Std(max)^2)^0.5*rand

    for i=1:length(ksi)
        if csy(1,i)>=csy(2,i)
            Xp(2,i)=ksi(2,i)*csy(2,i)+cay(2,i);
            Xp(1,i)=ksi(1,i)*(csy(1,i)^2-csy(2,i)^2)^(1/2)+(cay(1,i)-cay(2,i))+Xp(2,i);
        else
            Xp(1,i)=ksi(1,i)*csy(1,i)+cay(1,i);
            Xp(2,i)=ksi(2,i)*(csy(2,i)^2-csy(1,i)^2)^(1/2)-(cay(1,i)-cay(2,i))+Xp(1,i);
        end
    end

    % range control the Tmin, insure that Tmin is always small than Tmax
    for i=1:length(ksi)
        if Xp(1,i)<=Xp(2,i)
            Xp(2,i)=Xp(1,i)-abs(Xp(1,i))*0.2;
        end
    end
    Tmax=Xp(1,:);
    Tmin=Xp(2,:);
end

waitbar(0.9,h);
%
% store results in matrices of size [nYears 365]
%
gP=reshape(r,365,Y);
gP=gP';
clear r;
waitbar(0.93,h);

jj=find(gP>0);   % add precipitation threshold
gP(jj)=gP(jj)+PrecipThreshold;
clear jj;

gTmax=reshape(Tmax,365,Y);
gTmax=gTmax';
waitbar(0.96,h);
clear Tmax;
gTmin=reshape(Tmin,365,Y);
gTmin=gTmin';
clear Tmin;
%
% plot the first year results
figure
subplot(2,1,1)
bar(day,gP(1,:))
ylabel('daily precip. mm')
titre=['Generated daily precip of the first year'];
title(titre)
subplot(2,1,2)
plot(day,gTmax(1,:),day,gTmin(1,:),'r-')
ylabel('air temp,oC')
titre=['Generated daily Tmax and Tmin of the first year'];
title(titre)
legend('Tmax','Tmin','Location','Best')
%
% store results in files
%
save(filenameout,'gP','gTmax','gTmin');    % results stores in matrices [nYears 365]

waitbar(1,h);

close(h);