text:
function normalise = spm_cfg_norm
% SPM Configuration file for Spatial Normalisation
%__________________________________________________________________________
% Copyright (C) 2012-2016 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_cfg_norm.m 6952 2016-11-25 16:03:13Z guillaume $

%--------------------------------------------------------------------------
% biasreg Bias regularisation
%--------------------------------------------------------------------------
biasreg         = cfg_menu;
biasreg.tag     = 'biasreg';
biasreg.name    = 'Bias regularisation';
biasreg.help    = {
                   'MR images are usually corrupted by a smooth, spatially varying artifact that modulates the intensity of the image (bias). These artifacts, although not usually a problem for visual inspection, can impede automated processing of the images.'
                   ''
                   'An important issue relates to the distinction between intensity variations that arise because of bias artifact due to the physics of MR scanning, and those that arise due to different tissue properties.  The objective is to model the latter by different tissue classes, while modelling the former with a bias field. We know a priori that intensity variations due to MR physics tend to be spatially smooth, whereas those due to different tissue types tend to contain more high frequency information. A more accurate estimate of a bias field can be obtained by including prior knowledge about the distribution of the fields likely to be encountered by the correction algorithm. For example, if it is known that there is little or no intensity non-uniformity, then it would be wise to penalise large values for the intensity non-uniformity parameters. This regularisation can be placed within a Bayesian context, whereby the penalty incurred is the negative logarithm of a prior probability for any particular pattern of non-uniformity.'
                   'Knowing what works best should be a matter of empirical exploration.  For example, if your data has very little intensity non-uniformity artifact, then the bias regularisation should be increased.  This effectively tells the algorithm that there is very little bias in your data, so it does not try to model it.'
                   }';
biasreg.labels = {
                  'no regularisation (0)'
                  'extremely light regularisation (0.00001)'
                  'very light regularisation (0.0001)'
                  'light regularisation (0.001)'
                  'medium regularisation (0.01)'
                  'heavy regularisation (0.1)'
                  'very heavy regularisation (1)'
                  'extremely heavy regularisation (10)'
                  }';
biasreg.values = {
                  0
                  1e-05
                  0.0001
                  0.001
                  0.01
                  0.1
                  1
                  10
                  }';
biasreg.val    = {0.0001};

%--------------------------------------------------------------------------
% biasfwhm Bias FWHM
%--------------------------------------------------------------------------
biasfwhm        = cfg_menu;
biasfwhm.tag    = 'biasfwhm';
biasfwhm.name   = 'Bias FWHM';
biasfwhm.help   = {'FWHM of Gaussian smoothness of bias. If your intensity non-uniformity is very smooth, then choose a large FWHM. This will prevent the algorithm from trying to model out intensity variation due to different tissue types. The model for intensity non-uniformity is one of i.i.d. Gaussian noise that has been smoothed by some amount, before taking the exponential. Note also that smoother bias fields need fewer parameters to describe them. This means that the algorithm is faster for smoother intensity non-uniformities.'};
biasfwhm.labels = {
                   '30mm cutoff'
                   '40mm cutoff'
                   '50mm cutoff'
                   '60mm cutoff'
                   '70mm cutoff'
                   '80mm cutoff'
                   '90mm cutoff'
                   '100mm cutoff'
                   '110mm cutoff'
                   '120mm cutoff'
                   '130mm cutoff'
                   '140mm cutoff'
                   '150mm cutoff'
                   'No correction'
                   }';
biasfwhm.values = {
                   30
                   40
                   50
                   60
                   70
                   80
                   90
                   100
                   110
                   120
                   130
                   140
                   150
                   Inf
                   }';
biasfwhm.val    = {60};

%--------------------------------------------------------------------------
% write Save Bias Fields
%--------------------------------------------------------------------------
write        = cfg_menu;
write.tag    = 'write';
write.name   = 'Save Bias Fields';
write.help   = {'This is the option concerns whether to save the estimated bias fields. MR images are usually corrupted by a smooth, spatially varying artifact that modulates the intensity of the image (bias). These artifacts, although not usually a problem for visual inspection, can impede automated processing of the images.  The bias corrected version should have more uniform intensities within the different types of tissues.'};
write.labels = {
                'Save Nothing'
                'Save Bias Field'
                }';
write.values = {
                [0 0]
                [1 0]
                }';
write.val    = {[0 0]};

%--------------------------------------------------------------------------
% tpm Tissue probability map
%--------------------------------------------------------------------------
tpm         = cfg_files;
tpm.tag     = 'tpm';
tpm.name    = 'Tissue probability map';
tpm.help    = {
               'Select the tissue probability atlas. These should contain probability maps of all the various tissues found in the image data (such that probabilities are greater than or equal to zero, and they sum to one at each voxel. A nonlinear deformation field is estimated that best overlays the atlas on the individual subjects'' image.'
               }';
tpm.filter  = 'nifti';
tpm.ufilter = '.*';
tpm.num     = [1 1];
tpm.val     = {{fullfile(spm('dir'),'tpm','TPM.nii')}};
tpm.preview = @(f) spm_check_registration(char(f));

%--------------------------------------------------------------------------
% reg Warping Regularisation
%--------------------------------------------------------------------------
reg         = cfg_entry;
reg.tag     = 'reg';
reg.name    = 'Warping Regularisation';
reg.help    = {'The objective function for registering the tissue probability maps to the image to process, involves minimising the sum of two terms. One term gives a function of how probable the data is given the warping parameters. The other is a function of how probable the parameters are, and provides a penalty for unlikely deformations. Smoother deformations are deemed to be more probable. The amount of regularisation determines the tradeoff between the terms. Pick a value around one.  However, if your normalised images appear distorted, then it may be an idea to increase the amount of regularisation (by an order of magnitude). More regularisation gives smoother deformations, where the smoothness measure is determined by the bending energy of the deformations. '};
reg.strtype = 'r';
reg.num     = [1  5];
reg.val     = {[0 0.001 0.5 0.05 0.2]};
%%Eventually these values should be decreased (eg) to:
%reg.val     = {[0 0.0001 0.05 0.005 0.005]};

%--------------------------------------------------------------------------
% affreg Affine Regularisation
%--------------------------------------------------------------------------
affreg        = cfg_menu;
affreg.tag    = 'affreg';
affreg.name   = 'Affine Regularisation';
affreg.help   = {
                  'The procedure is a local optimisation, so it needs reasonable initial starting estimates. Images should be placed in approximate alignment using the Display function of SPM before beginning. A Mutual Information affine registration with the tissue probability maps (D''Agostino et al, 2004) is used to achieve approximate alignment. Note that this step does not include any model for intensity non-uniformity. This means that if the procedure is to be initialised with the affine registration, then the data should not be too corrupted with this artifact.If there is a lot of intensity non-uniformity, then manually position your image in order to achieve closer starting estimates, and turn off the affine registration.'
                  ''
                  'Affine registration into a standard space can be made more robust by regularisation (penalising excessive stretching or shrinking).  The best solutions can be obtained by knowing the approximate amount of stretching that is needed (e.g. ICBM templates are slightly bigger than typical brains, so greater zooms are likely to be needed). For example, if registering to an image in ICBM/MNI space, then choose this option.  If registering to a template that is close in size, then select the appropriate option for this.'
                  }';
affreg.labels = {
                 'No Affine Registration'
                 'ICBM space template - European brains'
                 'ICBM space template - East Asian brains'
                 'Average sized template'
                 'No regularisation'
                 }';
affreg.values = {
                 ''
                 'mni'
                 'eastern'
                 'subj'
                 'none'
                 }';
affreg.val    = {'mni'};

%--------------------------------------------------------------------------
% samp Sampling distance
%--------------------------------------------------------------------------
samp         = cfg_entry;
samp.tag     = 'samp';
samp.name    = 'Sampling distance';
samp.help    = {'This encodes the approximate distance between sampled points when estimating the model parameters. Smaller values use more of the data, but the procedure is slower and needs more memory. Determining the ``best'''' setting involves a compromise between speed and accuracy.'};
samp.strtype = 'r';
samp.num     = [1  1];
samp.val     = {3};

%--------------------------------------------------------------------------
% smo Smoothness
%--------------------------------------------------------------------------
smo         = cfg_entry;
smo.tag     = 'fwhm';
smo.name    = 'Smoothness';
smo.help    = {'For PET or SPECT, set this value to about 5 mm, or more if the images have smoother noise.  For MRI, you can usually use a value of 0 mm.  This is used to derive a fudge factor to account for correlations between neighbouring voxels.  Smoother data have more spatial correlations, rendering the assumptions of the model inaccurate.'};
smo.strtype = 'r';
smo.num     = [1  1];
smo.val     = {0};

%--------------------------------------------------------------------------
% write Deformation Fields
%--------------------------------------------------------------------------
write        = cfg_menu;
write.tag    = 'write';
write.name   = 'Deformation Fields';
write.help   = {'Deformation fields can be saved to disk, and used by the Deformations Utility. For spatially normalising images to MNI space, you will need the forward deformation, whereas for spatially normalising (eg) GIFTI surface files, you''ll need the inverse. It is also possible to transform data in MNI space on to the individual subject, which also requires the inverse transform. Deformations are saved as .nii files, which contain three volumes to encode the x, y and z coordinates.'};
write.labels = {
                'None'
                'Inverse'
                'Forward'
                'Inverse + Forward'
                }';
write.values = {
                [0 0]
                [1 0]
                [0 1]
                [1 1]
                }';
write.val    = {[0 0]};

%--------------------------------------------------------------------------
% eoptions Estimation Options
%--------------------------------------------------------------------------
eoptions      = cfg_branch;
eoptions.tag  = 'eoptions';
eoptions.name = 'Estimation Options';
eoptions.val  = {biasreg biasfwhm tpm affreg reg smo samp};
eoptions.help = {'Various settings for estimating deformations.'};

%--------------------------------------------------------------------------
% preserve Preserve
%--------------------------------------------------------------------------
preserve        = cfg_menu;
preserve.tag    = 'preserve';
preserve.name   = 'Preserve';
preserve.help   = {
                   'Preserve Concentrations: Spatially normalised images are not "modulated". The warped images preserve the intensities of the original images.'
                   ''
                   'Preserve Total: Spatially normalised images are "modulated" in order to preserve the total amount of signal in the images. Areas that are expanded during warping are correspondingly reduced in intensity.'
}';
preserve.labels = {
                   'Preserve Concentrations'
                   'Preserve Amount'
}';
preserve.values = {0 1};
preserve.def    = @(val)spm_get_defaults('normalise.write.preserve', val{:});

%--------------------------------------------------------------------------
% bb Bounding box
%--------------------------------------------------------------------------
bb         = cfg_entry;
bb.tag     = 'bb';
bb.name    = 'Bounding box';
bb.help    = {'The bounding box (in mm) of the volume which is to be written (relative to the anterior commissure).'};
bb.strtype = 'r';
bb.num     = [2 3];
bb.def     = @(val)spm_get_defaults('normalise.write.bb', val{:});

%--------------------------------------------------------------------------
% vox Voxel sizes
%--------------------------------------------------------------------------
vox         = cfg_entry;
vox.tag     = 'vox';
vox.name    = 'Voxel sizes';
vox.help    = {'The voxel sizes (x, y & z, in mm) of the written normalised images.'};
vox.strtype = 'r';
vox.num     = [1 3];
vox.def     = @(val)spm_get_defaults('normalise.write.vox', val{:});

%--------------------------------------------------------------------------
% interp Interpolation
%--------------------------------------------------------------------------
interp         = cfg_menu;
interp.tag     = 'interp';
interp.name    = 'Interpolation';
interp.help    = {
                  ['The method by which the images are sampled when ' ...
                  'being written in a different space. ' ...
                  '(Note that Inf or NaN values are treated as zero, ' ...
                  'rather than as missing data)']
                  '    Nearest Neighbour:'
                  '      - Fastest, but not normally recommended.'
                  '    Trilinear Interpolation:'
                  '      - OK for PET, realigned fMRI, or segmentations'
                  '    B-spline Interpolation:'
                  ['      - Better quality (but slower) interpolation' ...
                  '/* \cite{thevenaz00a}*/, especially with higher ' ...
                  'degree splines. Can produce values outside the ' ...
                  'original range (e.g. small negative values from an ' ...
                  'originally all positive image).']
}';
interp.labels = {
                 'Nearest neighbour'
                 'Trilinear'
                 '2nd Degree B-spline'
                 '3rd Degree B-Spline '
                 '4th Degree B-Spline '
                 '5th Degree B-Spline'
                 '6th Degree B-Spline'
                 '7th Degree B-Spline'
}';
interp.values = {0 1 2 3 4 5 6 7};
interp.def    = @(val)spm_get_defaults('normalise.write.interp', val{:});

%--------------------------------------------------------------------------
% wrap Wrapping
%--------------------------------------------------------------------------
wrap         = cfg_menu;
wrap.tag     = 'wrap';
wrap.name    = 'Wrapping';
wrap.help    = {
                'These are typically:'
                '    No wrapping: for PET or images that have already been spatially transformed. '
                '    Wrap in  Y: for (un-resliced) MRI where phase encoding is in the Y direction (voxel space).'
}';
wrap.labels  = {
               'No wrap'
               'Wrap X'
               'Wrap Y'
               'Wrap X & Y'
               'Wrap Z'
               'Wrap X & Z'
               'Wrap Y & Z'
               'Wrap X, Y & Z'
}';
wrap.values  = {[0 0 0] [1 0 0] [0 1 0] [1 1 0] [0 0 1] [1 0 1] [0 1 1]...
               [1 1 1]};
wrap.def     = @(val)spm_get_defaults('normalise.write.wrap', val{:});

%--------------------------------------------------------------------------
% prefix Filename Prefix
%--------------------------------------------------------------------------
prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the normalised image file(s). Default prefix is ''w''.'};
prefix.strtype = 's';
prefix.num     = [1 Inf];
prefix.def     = @(val)spm_get_defaults('normalise.write.prefix', val{:});

%--------------------------------------------------------------------------
% woptions Writing Options
%--------------------------------------------------------------------------
woptions      = cfg_branch;
woptions.tag  = 'woptions';
woptions.name = 'Writing Options';
woptions.val  = { bb vox interp prefix};
woptions.help = {'Various options for writing normalised images.'};

%--------------------------------------------------------------------------
% vol Image to Align
%--------------------------------------------------------------------------
vol         = cfg_files;
vol.tag     = 'vol';
vol.name    = 'Image to Align';
vol.help    = {
    'The image that the template (atlas) data is warped into alignment with.'
    'The result is a set of warps, which can be applied to this image, or any other image that is in register with it.'
    }';
vol.filter  = 'image';
vol.ufilter = '.*';
vol.num     = [1 1];
vol.preview = @(f) spm_check_registration(char(f));

%--------------------------------------------------------------------------
% def Parameter File
%--------------------------------------------------------------------------
def         = cfg_files;
def.tag     = 'def';
def.name    = 'Deformation Field';
def.help    = {[...
    'Deformations can be thought of as vector fields, and represented ',...
    'by three-volume images.  In SPM, deformation fields are saved in ',...
    'NIfTI format, with dimensions xdim x ydim x zdim x 1 x 3. ',...
    'Each voxel contains the x, y and z mm coordinates of where the deformation points.']};
def.filter  = 'nifti';
def.ufilter = 'y_.*\.nii$';
def.num     = [1 1];

%--------------------------------------------------------------------------
% resample Images to Write
%--------------------------------------------------------------------------
resample         = cfg_files;
resample.tag     = 'resample';
resample.name    = 'Images to Write';
resample.help    = {
    'These are the images for warping according to the estimated parameters.'
    'They can be any images that are in register with the image used to generate the deformation.'
    }';
resample.filter  = 'image';
resample.ufilter = '.*';
resample.num     = [1 Inf];
resample.preview = @(f) spm_check_registration(char(f));

%--------------------------------------------------------------------------
% subj Subject
%--------------------------------------------------------------------------
subj         = cfg_branch;
subj.tag     = 'subj';
subj.name    = 'Subject';
subj.val     = {vol};
subj.help    = {'Data for this subject. The same parameters are used within subject.'};

%--------------------------------------------------------------------------
% esubjs Data
%--------------------------------------------------------------------------
esubjs        = cfg_repeat;
esubjs.tag    = 'esubjs';
esubjs.name   = 'Data';
esubjs.help   = {'List of subjects. Images of each subject should be warped differently.'};
esubjs.values = {subj};
esubjs.num    = [1 Inf];

%--------------------------------------------------------------------------
% subj Subject
%--------------------------------------------------------------------------
subj      = cfg_branch;
subj.tag  = 'subj';
subj.name = 'Subject';
subj.val  = {def resample};
subj.help = {'Data for this subject. The same parameters are used within subject.'};

%--------------------------------------------------------------------------
% wsubjs Data
%--------------------------------------------------------------------------
wsubjs        = cfg_repeat;
wsubjs.tag    = 'wsubjs';
wsubjs.name   = 'Data';
wsubjs.help   = {'List of subjects. Images of each subject should be warped differently.'};
wsubjs.values = {subj};
wsubjs.num    = [1 Inf];

%--------------------------------------------------------------------------
% subj Subject
%--------------------------------------------------------------------------
subj      = cfg_branch;
subj.tag  = 'subj';
subj.name = 'Subject';
subj.val  = {vol resample};
subj.help = {'Data for this subject. The same parameters are used within subject.'};

%--------------------------------------------------------------------------
% ewsubjs Data
%--------------------------------------------------------------------------
ewsubjs        = cfg_repeat;
ewsubjs.tag    = 'ewsubjs';
ewsubjs.name   = 'Data';
ewsubjs.help   = {'List of subjects. Images of each subject should be warped differently.'};
ewsubjs.values = {subj};
ewsubjs.num    = [1 Inf];

%--------------------------------------------------------------------------
% est Segment
%--------------------------------------------------------------------------
est      = cfg_exbranch;
est.tag  = 'est';
est.name = 'Normalise: Estimate';
est.val  = {esubjs eoptions};
est.help = {
    'Spatial normalisation performed via the segmentation routine.'
    ''
    'The algorithm (which was known as ``New Segment'''' in SPM8) is essentially the same as that described in the Unified Segmentation paper /* \cite{ashburner05}*/, except for (i) a slightly different treatment of the mixing proportions, (ii) the use of an improved registration model, (iii) the ability to use multi-spectral data, (iv) an extended set of tissue probability maps, which allows a different treatment of voxels outside the brain.'
    ''
    'If you encounter problems with spatial normalisation, it is advisable to use the Check reg button to see how well aligned the original data are with the MNI-space templates released with SPM.  If mis-alignment is greater than about 3cm and 15 degrees, you could try to manually re-position the images prior to attempting to align them.  This may be done using the Display button.'
    }';
est.prog = @spm_run_norm;
est.vout = @vout_est;

%--------------------------------------------------------------------------
% write Normalise: Write
%--------------------------------------------------------------------------
write      = cfg_exbranch;
write.tag  = 'write';
write.name = 'Normalise: Write';
write.val  = {wsubjs woptions};
write.help = {
    'Apply previously estimated warps (stored in ``y_''''imagename``_sn.mat'''' files) to series of images.'};
write.prog = @spm_run_norm;
write.vout = @vout_write;

%--------------------------------------------------------------------------
% estwrite Normalise: Estimate & Write
%--------------------------------------------------------------------------
estwrite      = cfg_exbranch;
estwrite.tag  = 'estwrite';
estwrite.name = 'Normalise: Estimate & Write';
estwrite.val  = {ewsubjs eoptions woptions};
estwrite.help = {
    'Compute the warp that best aligns the template (atlas) to the individual''s image, invert it and write the result to the file `y_''imagename''.nii''.'
    'This option also allows the contents of the `y_''imagename''.nii'' files to be applied to a series of images.'
    ''
    'Note that if you encounter problems with spatial normalisation, it is often advisable to use the Check reg button to see how well aligned the original data are with the MNI-space templates released with SPM.  If mis-alignment is greater than about 3cm and 15 degrees, you could try to manually re-position the images.  This may be done using the Display button.'
    };
estwrite.prog = @spm_run_norm;
estwrite.vout = @vout_estwrite;

%--------------------------------------------------------------------------
% normalise Normalise
%--------------------------------------------------------------------------
normalise        = cfg_choice;
normalise.tag    = 'normalise';
normalise.name   = 'Normalise';
normalise.help   = {...
    'There are two components to spatial normalisation: There is the estimation part, ',...
    'whereby a deformation is estimated by deforming template data to match an ',...
    'individual scan; And there is the actual writing of the spatially normalised ',...
    'images, using the previously estimated deformation.',...
    'This is a vanilla approach to spatial normalisation.  ',...
    'It is not generally recommended for morphometric studies, or other studies of ',...
    'differences among populations. ',...
    'The reason is that the atlas data will differ systematically from the data under study, ',...
    'which is likely to lead to an inherently biased set of findings.' };
normalise.values = {est write estwrite};

%==========================================================================
function dep = vout_est(job)
for k=1:numel(job.subj)
    dep(k)            = cfg_dep;
    dep(k).sname      = sprintf('Deformation (Subj %d)',k);
    dep(k).src_output = substruct('()',{k},'.','def');
    dep(k).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
end

%==========================================================================
function dep = vout_write(job)
for k=1:numel(job.subj)
    dep(k)            = cfg_dep;
    dep(k).sname      = sprintf('Normalised Images (Subj %d)',k);
    dep(k).src_output = substruct('()',{k},'.','files');
    dep(k).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
end

%==========================================================================
function dep = vout_estwrite(job)
depe = vout_est(job);
depw = vout_write(job);
dep = [depe depw];

author: spm


text:
function img = gen_img_mesh(size, res)
    img = zeros(size);
    qs = 1 : res : size(1);
    for i = 1 : numel(qs)
%         qs(i)
        img(ceil(qs(i)), :) = 1;
    end
    qs = 1 : res : size(2);
    for i = 1 : numel(qs)
        img(:, ceil(qs(i))) = 1;
    end
    img(end, :) = 1;
    img(:, end) = 1;
end

repo: pTVreg


text:
function DI = detectChange(obj, t1, t2)
DI = abs(double(t1) - double(t2));
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
function [tfr,t,f] = tfrpmh(x,t,N,h,trace);
%TFRPMH  Pseudo Margenau-Hill time-frequency distribution.
%       [TFR,T,F]=TFRPMH(X,T,N,H,TRACE) computes the Pseudo Margenau-Hill
%       distribution of a discrete-time signal X, or the
%       cross Pseudo Margenau-Hill representation between two signals.
%
%       X     : signal if auto-PMH, or [X1,X2] if cross-PMH.
%       T     : time instant(s)          (default : 1:length(X)).
%       N     : number of frequency bins (default : length(X)).
%       H     : frequency smoothing window, H(0) being forced to 1
%                                        (default : Hamming(N/4)).
%       TRACE : if nonzero, the progression of the algorithm is shown
%                                        (default : 0).
%       TFR   : time-frequency representation. When called without
%               output arguments, TFRPMH runs TFRQVIEW.
%       F     : vector of normalized frequencies.
%
%       Example :
%        sig=fmlin(128,0.1,0.4); t=1:128;
%        h=tftb_window(63,'Kaiser'); tfrpmh(sig,t,128,h,1);
%
%       See also all the time-frequency representations listed in
%        the file CONTENTS (TFR*)

%       F. Auger, May-August 1994, July 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

[xrow,xcol] = size(x);
if (nargin < 1),
 error('At least 1 parameter is required');
elseif nargin<=2,
 N=xrow;
end;

hlength=floor(N/4);
hlength=hlength+1-rem(hlength,2);

if (nargin == 1),
 t=1:xrow; h = tftb_window(hlength); trace=0;
elseif (nargin == 2 | nargin == 3),
 h = tftb_window(hlength); trace = 0;
elseif (nargin == 4),
 trace = 0;
end;

if (N<0),
 error('N must be greater than zero');
end;
[trow,tcol] = size(t);
if (xcol==0)|(xcol>2),
 error('X must have one or two columns');
elseif (trow~=1),
 error('T must only have one row');
elseif (2^nextpow2(N)~=N), %(rem(log(N)/log(2),1)~=0),
 fprintf('For a faster computation, N should be a power of two\n');
end;

[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1);
if (hcol~=1)|(rem(hrow,2)==0),
 error('H must be a smoothing window with odd length');
end;

tfr= zeros (N,tcol) ;
if trace, disp('Pseudo Margenau-Hill distribution'); end;
for icol=1:tcol,
 ti= t(icol); tau=-min([round(N/2)-1,Lh,xrow-ti]):min([round(N/2)-1,Lh,ti-1]);
 indices= rem(N+tau,N)+1;
 if trace, disprog(icol,tcol,10); end;
 tfr(indices,icol)=h(Lh+1+tau).*x(ti,1).*conj(x(ti-tau,xcol));
end;
if trace, fprintf('\n'); end;
tfr= real(fft(tfr));

if (nargout==0),
 tfrqview(tfr,x,t,'tfrpmh',h);
elseif (nargout==3),
 if rem(N,2)==0,
  f=[0:N/2-1 -N/2:-1]'/N;
 else
  f=[0:(N-1)/2 -(N-1)/2:-1]'/N;
 end;
end;

author: HeLiangHIT


text:
function X = tenzeros(varargin)
%TENZEROS Create zeros tensor.
%
%   X = TENZEROS(SZ) forms a tensor of size SZ with all zeros.
%
%   TENZEROS(SZ) is equivalent to TENSOR(ZEROS(SZ(1),SZ(2),...),SZ).
%
%   See also TENSOR, ZEROS.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

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

if nargin == 1
    sz = varargin{1};
else
    sz = cell2mat(varargin);
end

if isempty(sz)
    X = tensor;
    return;
end

if nargin == 2
    order = sz;
    dim = varargin{1};
    sz = dim * ones(1,order);
end

data = zeros([sz 1 1]);
X = tensor(data,sz);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function layer = genNetworkLSTMEncoderDecoderScale2(para)

layer = genNetworkLSTMEncoderDecoder2(para);

layerScale.name = 'Input';
layerScale.inputIdx = 3;
layerScale.dim = [1 1];
layerScale.next = 1;

layer = [layer(1:end-1) layerScale layer(end)];

layer{end-3}.next = 3;
layer{end-2}.next = 2;
layer{end}.prev = [-3 -2 -1];
% layer{end} = rmfield(layer{end}, 'costFrameSelection');

end

author: singaxiong

repo: SignalGraph


text:
function b = r8ge_ml ( n, a_lu, pivot, x, job )

%*****************************************************************************80
%
%% R8GE_ML computes A * x or A' * x, using R8GE_FA factors.
%
%  Discussion:
%
%    The R8GE storage format is used for a general M by N matrix.  A storage
%    space is made for each logical entry.  The two dimensional logical
%    array is mapped to a vector, in which storage is by columns.
%
%    It is assumed that R8GE_FA has overwritten the original matrix
%    information by LU factors.  R8GE_ML is able to reconstruct the
%    original matrix from the LU factor data.
%
%    R8GE_ML allows the user to check that the solution of a linear
%    system is correct, without having to save an unfactored copy
%    of the matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, real A_LU(N,N), the LU factors from R8GE_FA.
%
%    Input, integer PIVOT(N), the pivot vector computed by R8GE_FA.
%
%    Input, real X(N), the vector to be multiplied.
%
%    Input, integer JOB, specifies the operation to be done:
%    JOB = 0, compute A * x.
%    JOB nonzero, compute A' * X.
%
%    Output, real B(N), the result of the multiplication.
%
  b(1:n) = x(1:n);

  if ( job == 0 )
%
%  Y = U * X.
%
    for j = 1 : n
      b(1:j-1) = b(1:j-1) + a_lu(1:j-1,j)' * b(j);
      b(j) = a_lu(j,j) * b(j);
    end
%
%  B = PL * Y = PL * U * X = A * x.
%
    for j = n-1 : -1 : 1

      b(j+1:n) = b(j+1:n) - a_lu(j+1:n,j)' * b(j);
      k = pivot(j);

      if ( k ~= j )
        t    = b(k);
        b(k) = b(j);
        b(j) = t;
      end

    end

  else
%
%  Y = (PL)' * X:
%
    for j = 1 : n-1

      k = pivot(j);

      if ( k ~= j )
        t    = b(k);
        b(k) = b(j);
        b(j) = t;
      end

      b(j) = b(j) - sum ( b(j+1:n) * a_lu(j+1:n,j) );

    end
%
%  B = U' * Y = ( PL * U )' * X = A' * X.
%
    for i = n : -1 : 1
      b(i+1:n) = b(i+1:n) + b(i) * a_lu(i,i+1:n);
      b(i) = b(i) * a_lu(i,i);
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%
% Set the default parameters for each model
%
function AMICO_SetModel( model )

    global CONFIG

    CONFIG.model   = [];

    CONFIG.OPTIMIZATION = [];
    CONFIG.OPTIMIZATION.LS_param    = optimset('TolX',1e-4);
    CONFIG.OPTIMIZATION.SPAMS_param = [];

    % Call the specific model constructor
    modelClass = str2func( ['AMICO_' upper(model)] );
    if exist([func2str(modelClass) '.m'],'file')
        CONFIG.model = modelClass();
    else
        error( '[AMICO_SetModel] Model "%s" not recognized', model )
    end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function r8vec_sort_heap_a_test ( )

%*****************************************************************************80
%
%% R8VEC_SORT_HEAP_A_TEST tests R8VEC_SORT_HEAP_A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 20;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_SORT_HEAP_A_TEST\n' );
  fprintf ( 1, '  R8VEC_SORT_HEAP_A ascending sorts an R8VEC.\n' );

  b = 0.0;
  c = 3.0 * n;
  seed = 123456789;

  [ a, seed ] = r8vec_uniform_ab ( n, b, c, seed );

  r8vec_print_some ( n, a, 1, 10, '  Original array:' );

  a = r8vec_sort_heap_a ( n, a );

  r8vec_print_some ( n, a, 1, 10, '  Ascending sorted array:' );

  return
end


text:
function [colorCodes, signChange, errorFlag]= getColorCodes(referenceFluxes, fluxes, maxChange, defaultColor, upColor, downColor)
% getColorCodes
%	Calculates the coloring for a number of fluxes by comparing them to
%   reference fluxes.
%
%   referenceFluxes     vector of reference fluxes
%   fluxes              vector of fluxes. The number of elements in fluxes
%                       and referenceFluxes must be equal
%   maxChange           the logfold increase or decrease that corresponds
%                       to full negative or full positive coloring. Must
%                       be a positive value (opt, default 1)
%   defaultColor        a color in Matlab format to be used if there are no
%                       changes between the fluxes. This color is also used to
%                       calculate the transition between the colors for up and
%                       down regulated fluxes (opt, default [1 1 1])
%   upColor             a color in Matlab format to be used if the flux is
%                       larger than the reference flux (opt, default [0 1 0])
%   downColor           a color in Matlab format to be used if the flux is
%                       smaller than the reference flux (opt, default [1 0 0])

%   colorCodes          array list of colors in Matlab format in the same
%                       order as the fluxes
%   signChange          array list of boolean values where true indicates
%                       that there has been a sign change between the
%                       fluxes. Reactions with sign changes are not
%                       colored, but rather marked in another way. The
%                       order correspondes to the order of the fluxes
%   errorFlag           true if there has been an error
%
%   Usage: [colorCodes, signChange, errorFlag]=getColorCodes(referenceFluxes,...
%           fluxes, maxChange, defaultColor, upColor, downColor)

if nargin<6
    downColor=[1 0 0];
end
if nargin<5
    upColor=[0 1 0];
end
if nargin<4
    defaultColor=[1 1 1];
end
if nargin<3
    maxChange=1;
end

%Checks that the flux vector and the reference flux vector have the same
%number of elements
if length(fluxes)~=length(referenceFluxes)
    colorCodes={};
    signChange={};
    errorFlag=1;
    fprintf('fluxes and referenceFluxes must have the same dimensions');
    return;
end

%Loops through the fluxes. If the reference flux is 0, then mark as
%upColor, if the flux is 0 then mark as downColor, and if both fluxes are 0
%then mark as defaultColor
for i=1:length(fluxes)
    signChange{i}=false;

    if referenceFluxes(i)==0 || fluxes(i)==0
        if referenceFluxes(i)==0 && fluxes(i)==0
            colorCodes{i}=defaultColor;
        else
            if referenceFluxes(i)==0
                colorCodes{i}=upColor;
            else
                colorCodes{i}=downColor;
            end
        end
    else
        %At the moment a negative flux that gets more negative is counted
        %as being upregulated. This might be counter intuitive.
        logvalue=log10(abs(fluxes(i))/abs(referenceFluxes(i)));

        %If there is a sign change
        if fluxes(i)*referenceFluxes(i)<0
            colorCodes{i}=defaultColor;
            signChange{i}=true;
        else
            colorCodes{i}=getColor(logvalue, defaultColor, upColor, downColor, maxChange);
        end
    end
end
end

function colorValue=getColor(logvalue, defaultColor, upColor, downColor, maxChange)
% getColor
%   Calculates the color for a specified logvalue.
%
%   logvalue            the logfold increase or desrease of a flux
%   compared to the
%                       corresponding reference flux. Must be a
%                       positive value
%   defaultColor        a color in Matlab format to be used if
%   there are no
%                       changes between the fluxes. This color is
%                       also used to calculate the transition
%                       between the colors for up and down
%                       regulated fluxes
%   upColor             a color in Matlab format to be used if the
%   flux is
%                       larger than the reference flux
%   downColor           a color in Matlab format to be used if the
%   flux is
%                       smaller than the reference flux
%   maxChange           the logfold increase or decrease that
%   corresponds
%                       to full negative or full positive coloring
%
%   colorValue          vector with the calculated color
%
%   Usage: colorValue=getColor(logvalue, defaultColor, upColor,
%   downColor, maxChange)

%If the flux has decreased
if logvalue<0
    %If the flux is lower than 10^-maxChange of the original then
    %color as downColor
    if logvalue<-maxChange
        colorValue=downColor;
    else
        %The color is linear from defaultColor to downColor
        colorValue=[defaultColor(1)+(downColor(1)-defaultColor(1))*logvalue/(-maxChange)...
            defaultColor(2)+(downColor(2)-defaultColor(2))*logvalue/(-maxChange)...
            defaultColor(3)+(downColor(3)-defaultColor(3))*logvalue/(-maxChange)];
    end
    %If it has increased
else
    %If the flux is higher than 10^maxChange times the original
    %then color green
    if logvalue>maxChange
        colorValue=upColor;
    else
        %The color is linear from defaultColor to upColor
        colorValue=[defaultColor(1)+(upColor(1)-defaultColor(1))*logvalue/(maxChange)...
            defaultColor(2)+(upColor(2)-defaultColor(2))*logvalue/(maxChange)...
            defaultColor(3)+(upColor(3)-defaultColor(3))*logvalue/(maxChange)];
    end
end
end

author: SysBioChalmers

repo: RAVEN


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%|             francois.alouges@polytechnique.edu                         |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtDomEdge.m                                  |
%|    #    |   VERSION    : 0.50                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 25.11.2018                                    |
%|  / 0 \  |   LAST MODIF : 25.11.2018                                    |
%| ( === ) |   SYNOPSIS   : Triangular quadrature and integration         |
%|  `---'  |                                                              |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Parameters
Nvtx = 1e3;
rho  = 1;
tol  = 1e-3;

% Square mesh
mesh = mshSegment(Nvtx,rho);

% Domain
omega = dom(mesh,3);

% Finite element space
u = fem(mesh,'P1');
v = fem(mesh,'P1');

% Graphical representation
figure
plot(mesh)
hold on
plotNrm(mesh)
plot(omega)
plot(u,'or')
axis equal
alpha(0.1)
view(0,90)

% Numerical functions
Fx   = @(X) ones(size(X,1),1);
Fx3  = {Fx,Fx,Fx};
Fxy  = @(X,Y) (X(:,1)==Y(:,1));
Fxy3 = {Fxy,Fxy,Fxy};
Gxy  = @(X,Y) 1./(4*pi) * femGreenKernel(X,Y,'[exp(ikr)/r]',5);
Gxy3 = {Gxy,Gxy,Gxy};

%%%%%%%%%%%%%%% 2 ARGUMENTS %%%%%%%%%%%%%%%
disp('=========== 2 ARGUMENTS  ============')

% \int_{mesh(x)} f(x) dx
ref = integral(omega,Fx);
abs(ref-1)

% \int_{mesh(x)} f3(x) dx
sol = integral(omega,Fx3);
norm(sol{1}-ref,'inf')

%%%%%%%%%%%%%%% 3 ARGUMENTS %%%%%%%%%%%%%%%
disp('=========== 3 ARGUMENTS  ============')

% \int_{mesh(y)} f(x,y) dy
ref = integral(omega.qud,omega,Fxy);
abs(sum(ref,1)-1)

% \int_{mesh(y)} f3(x,y) dy
sol = integral(omega.qud,omega,Fxy3);
norm(sol{1}-ref,'inf')

%-------------------------------------

% \int_{mesh(x)} f(x,y) dx
ref = integral(omega,omega.qud,Fxy);
abs(sum(ref,2)-1)

% \int_{mesh(y)} f3(x,y) dx
sol = integral(omega,omega.qud,Fxy3);
norm(sol{1}-ref,'inf')

%-------------------------------------

% \int_{mesh(x)} f(x) psi(x) dx
ref = integral(omega,Fx,v);
abs(sum(ref,2) - 1)

% \int_{mesh(x)} f3(x) psi(x) dx
sol = integral(omega,Fx3,v);
norm(sol{1}-ref,'inf')

% \int_{mesh(x)} f(x) ntimes(psi(x)) dx
sol = integral(omega,Fx,ntimes(v));
norm(sol{2}-ref,'inf')

% \int_{mesh(x)} f3(x) ntimes(psi(x)) dx
sol = integral(omega,Fx3,ntimes(v));
norm(sol-ref,'inf')

%-------------------------------------

% \int_{mesh(x)} psi(x)' f(x)  dx
ref = integral(omega,u,Fx);
abs(sum(ref,1) - 1)

% \int_{mesh(x)} psi(x)' f3(x)  dx
sol = integral(omega,u,Fx3);
norm(sol{1}-ref,'inf')

% \int_{mesh(x)} ntimes(psi(x))' f(x)  dx
sol = integral(omega,ntimes(u),Fx);
norm(sol{2}-ref,'inf')

% \int_{mesh(x)} ntimes(psi(x))' f3(x)  dx
sol = integral(omega,ntimes(u),Fx3);
norm(sol-ref,'inf')

%-------------------------------------

% \int_{mesh(x)} psi(x)' psi(x) dx
ref = integral(omega,u,v);
abs(sum(sum(ref,1),2) - 1)

% \int_{mesh(x)} ntimes(psi(x))' psi(x) dx
sol = integral(omega,ntimes(u),v);
norm(sol{2}-ref,'inf')

% \int_{mesh(x)} psi(x)' ntimes(psi(x)) dx
sol = integral(omega,u,ntimes(v));
norm(sol{2}-ref,'inf')

% \int_{mesh(x)} ntimes(psi(x))' ntimes(psi(x)) dx
sol = integral(omega,ntimes(u),ntimes(v));
norm(sol-ref,'inf')

%%%%%%%%%%%%%%% 4 ARGUMENTS %%%%%%%%%%%%%%%
disp('=========== 4 ARGUMENTS  ============')

% \int_{mesh(x)} psi(x)' f(x) psi(x) dx
ref = integral(omega,u,Fx,v);
abs(sum(sum(ref,1),2) - 1)

% \int_{mesh(x)} ntimes(psi(x))' f(x) psi(x) dx
sol = integral(omega,ntimes(u),Fx,v);
norm(sol{2}-ref,'inf')

% \int_{mesh(x)} psi(x)' f3(x) psi(x) dx
sol = integral(omega,u,Fx3,v);
norm(sol{2}-ref,'inf')

% \int_{mesh(x)} psi(x)' f(x) ntimes(psi(x)) dx
sol = integral(omega,u,Fx,ntimes(v));
norm(sol{2}-ref,'inf')

% \int_{mesh(x)} ntimes(psi(x))' f3(x) psi(x) dx
sol = integral(omega,ntimes(u),Fx3,v);
norm(sol-ref,'inf')

% \int_{mesh(x)} ntimes(psi(x))' f(x) ntimes(psi(x)) dx
sol = integral(omega,ntimes(u),Fx,ntimes(v));
norm(sol-ref,'inf')

% \int_{mesh(x)} (psi(x))' f3(x) ntimes(psi(x)) dx
sol = integral(omega,u,Fx3,ntimes(v));
norm(sol-ref,'inf')

% \int_{mesh(x)} ntimes(psi(x))' (f3(x)) ntimes(psi(x)) dx
sol = integral(omega,ntimes(u),Fx3,ntimes(v));
norm(sol,'inf')

%-------------------------------------

% \int_{mesh(y)} f(x,y) psi(y) dy
ref = integral(omega.qud,omega,Fxy,v);
abs(sum(sum(ref,1),2) - 1)

% \int_{mesh(y)} f(x,y) ntimes(psi(y)) dy
sol = integral(omega.qud,omega,Fxy,ntimes(v));
norm(sol{2}-ref,'inf')

% \int_{mesh(y)} f(3x,y) psi(y) dy
sol = integral(omega.qud,omega,Fxy3,v);
norm(sol{2}-ref,'inf')

% \int_{mesh(y)} f3(x,y) ntimes(psi(y)) dy
sol = integral(omega.qud,omega,Fxy3,ntimes(v));
norm(sol-ref,'inf')

%-------------------------------------

% \int_{mesh(x)} psi(x)' f(x,y) dx
ref = integral(omega,omega.qud,u,Fxy);
abs(sum(sum(ref,1),2) - 1)

% \int_{mesh(x)} psi(x)' f3(x,y) dx
sol = integral(omega,omega.qud,u,Fxy3);
norm(sol{2}-ref,'inf')

% \int_{mesh(x)} ntimes(psi(x))' f(x,y) dx
sol = integral(omega,omega.qud,ntimes(u),Fxy);
norm(sol{2}-ref,'inf')

% \int_{mesh(x)} ntimes(psi(x))' f3(x,y) dx
sol = integral(omega,omega.qud,ntimes(u),Fxy3);
norm(sol-ref,'inf')

%%%%%%%%%%%%%% 5 ARGUMENTS %%%%%%%%%%%%%%%
disp('=========== 5 ARGUMENTS  ============')

% \int_{mesh(y)} G(x,y) psi(y) dy
ref = integral(mesh.ctr,omega,Gxy,v);
sol = integral(mesh.ctr,omega,Gxy,v,tol);
norm(full(sol)-ref,'inf')./norm(ref,'inf')

% % \int_{mesh(y)} G3(x,y) psi(y) dy
sol = integral(mesh.ctr,omega,Gxy3,v,tol);
norm(full(sol{1})-ref,'inf')./norm(ref,'inf')

% \int_{mesh(y)} G(x,y) ntimes(psi(y)) dy
sol = integral(mesh.ctr,omega,Gxy,ntimes(v),tol);
norm(full(sol{2})-ref,'inf')./norm(ref,'inf')

% \int_{mesh(y)} G3(x,y) ntimes(psi(y)) dy
sol = integral(mesh.ctr,omega,Gxy3,ntimes(v),tol);
norm(full(sol)-ref,'inf')./norm(ref,'inf')

%-------------------------------------

% \int_{mesh(x)} psi(x)' G(x,y) dx
ref = integral(omega,mesh.ctr,u,Gxy);
sol = integral(omega,mesh.ctr,u,Gxy,tol);
norm(full(sol)-ref,'inf')./norm(ref,'inf')

% \int_{mesh(x)} psi(x)' G3(x,y) dx
sol = integral(omega,mesh.ctr,u,Gxy3,tol);
norm(full(sol{1})-ref,'inf')./norm(ref,'inf')

% \int_{mesh(x)} ntimes(psi(x))' G(x,y) dx
sol = integral(omega,mesh.ctr,ntimes(u),Gxy,tol);
norm(full(sol{2})-ref,'inf')./norm(ref,'inf')

% \int_{mesh(x)} ntimes(psi(x))' G3(x,y) dx
sol = integral(omega,mesh.ctr,ntimes(u),Gxy3,tol);
norm(full(sol)-ref,'inf')./norm(ref,'inf')

%%%%%%%%%%%%%%% 6 ARGUMENTS %%%%%%%%%%%%%%%
disp('=========== 6 ARGUMENTS  ============')

% \int_{mesh(x)} \int_{mesh(y)} psi(x)' G(x,y) psi(y) dx dy
ref = integral(omega,omega,u,Gxy,v);
sol = integral(omega,omega,u,Gxy,v,tol);
norm(full(sol)-ref,'inf')./norm(ref,'inf')

% \int_{mesh(x)} \int_{mesh(y)} ntimes(psi(x))' G(x,y) psi(y) dx dy
ref = integral(omega,omega,ntimes(u),Gxy,v);
sol = integral(omega,omega,ntimes(u),Gxy,v,tol);
norm(full(sol{2})-ref{2},'inf')./norm(ref{2},'inf')

% \int_{mesh(x)} \int_{mesh(y)} psi(x)' G3(x,y) psi(y) dx dy
ref = integral(omega,omega,u,Gxy3,v);
sol = integral(omega,omega,u,Gxy3,v,tol);
norm(full(sol{2})-ref{2},'inf')./norm(ref{2},'inf')

% \int_{mesh(x)} \int_{mesh(y)} psi(x)' G(x,y) ntimes(psi(y)) dx dy
ref = integral(omega,omega,u,Gxy,ntimes(v));
sol = integral(omega,omega,u,Gxy,ntimes(v),tol);
norm(full(sol{2})-ref{2},'inf')./norm(ref{2},'inf')

% \int_{mesh(x)} \int_{mesh(y)} ntimes(psi(x))' G3(x,y) psi(y) dx dy
ref = integral(omega,omega,ntimes(u),Gxy3,v);
sol = integral(omega,omega,ntimes(u),Gxy3,v,tol);
norm(full(sol)-ref,'inf')./norm(ref,'inf')

% \int_{mesh(x)} \int_{mesh(y)} psi(x)' G3(x,y) ntimes(psi(y)) dx dy
ref = integral(omega,omega,u,Gxy3,ntimes(v));
sol = integral(omega,omega,u,Gxy3,ntimes(v),tol);
norm(full(sol)-ref,'inf')./norm(ref,'inf')

% \int_{mesh(x)} \int_{mesh(y)} ntimes(psi(x))' G(x,y) ntimes(psi(y)) dx dy
ref = integral(omega,omega,ntimes(u),Gxy,ntimes(v));
sol = integral(omega,omega,ntimes(u),Gxy,ntimes(v),tol);
norm(full(sol)-ref,'inf')./norm(ref,'inf')

% \int_{mesh(x)} \int_{mesh(y)} ntimes(psi(x))' G3(x,y) ntimes(psi(y)) dx dy
ref = integral(omega,omega,ntimes(u),Gxy3,ntimes(v));
sol = integral(omega,omega,ntimes(u),Gxy3,ntimes(v),tol);
norm(full(sol)-ref,'inf')

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% Extract critical features (junctions and endpoints).
%
% Rougly follows procedure of Liu et al. "Identification of Fork Points..." IEEE
% TAPMI
%
% For a circular stroke, there are no "features" by this definition.
% Thus, for any partition of the pixels into regions, there should
% be at least one feature for the tracing algorithm to find.
%
% Input
%  T: [n x n boolean] thinned image.
%    images are binary, where true means "black"
%
% Output
%  SN: [n x n boolean] extracted features.
function SN = extract_junctions(T)

    SE = bwmorph(T,'endpoints');
    SB = T; % black pixels
    sz = size(T,1);

    lutS3 = makelut( @(P)fS3(P) , 3);
    S3 = applylut(T,lutS3);

    % final criteria
    SN = SE | (SB & S3);

    % Check to see that each connected component has a feature.
    % This is necessary to process circles in the image.
    CC = bwconncomp(T,8);
    nCC = CC.NumObjects;
    for c=1:nCC

       pid = CC.PixelIdxList{c};

       % We have a circle. Circles are generally drawn from the
       % top, we choose the top pixel here
       if sum(SN(pid))==0
          [irow,icol] = ind2sub(sz,pid);
          sel = argmin(irow);
          SN(pid(sel)) = true;
       end

    end

end

% See Liu et al.
function Y=fS3(P)
    sz = size(P);
    assert(isequal(sz,[3 3]));

    % Get cross number
    NC = fNC(P);

    % Count black pixels
    PM = P;
    PM(2,2) = false;
    NB = sum(PM(:));

    % Criteria
    Y = (NC >= 3-eps) || (NB >= 4-eps);
end

% See Liu et al.
function Y=fNC(P)
    sum = 0;
    for i=0:7
        sum = sum + abs( P(fIP(i+1)) - P(fIP(i)) );
    end
    Y = sum./2;
end

% See Liu et al.
function newlindx = fIP(lindx)
    switch lindx
        case {0,8}
            i=1; j=2;
        case 1
            i=1; j=3;
        case 2
            i=2; j=3;
        case 3
            i=3; j=3;
        case 4
            i=3; j=2;
        case 5
            i=3; j=1;
        case 6
            i=2; j=1;
        case 7
            i=1; j=1;
    end
    newlindx = sub2ind([3 3],i,j);
end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function x = r83_cr_sl ( n, a_cr, b )

%*****************************************************************************80
%
%% R83_CR_SL solves a real linear system factored by R83_CR_FA.
%
%  Discussion:
%
%    The matrix A must be tridiagonal.  R83_CR_FA is called to compute the
%    LU factors of A.  It does so using a form of cyclic reduction.  If
%    the factors computed by R83_CR_FA are passed to R83_CR_SL, then one or many
%    linear systems involving the matrix A may be solved.
%
%    Note that R83_CR_FA does not perform pivoting, and so the solution
%    produced by R83_CR_SL may be less accurate than a solution produced
%    by a standard Gauss algorithm.  However, such problems can be
%    guaranteed not to occur if the matrix A is strictly diagonally
%    dominant, that is, if the absolute value of the diagonal coefficient
%    is greater than the sum of the absolute values of the two off diagonal
%    coefficients, for each row of the matrix.
%
%  Example:
%
%    Here is how a R83 matrix of order 5 would be stored:
%
%       *  A12 A23 A34 A45
%      A11 A22 A33 A44 A55
%      A21 A32 A43 A54  *
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 March 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Roger Hockney,
%    A fast direct solution of Poisson's equation using Fourier Analysis,
%    Journal of the ACM,
%    Volume 12, Number 1, pages 95-113, January 1965.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, real A_CR(3,2*N+1), factorization information computed by R83_CR_FA.
%
%    Input, real B(N), the right hand side vector.
%
%    Output, real X(N), the solution of the linear system.
%
  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R83_CR_SL - Fatal error!\n' );
    fprintf ( 1, '  Nonpositive N = %d\n', n );
    return
  end

  if ( n == 1 )
    x(1) = a_cr(2,2) * b(1);
    return
  end
%
%  Set up RHS.
%
  rhs(1) = 0.0;
  rhs(2:n+1) = b(1:n);
  rhs(n+2:2*n+1) = 0.0;

  il = n;
  ndiv = 1;
  ipntp = 0;

  while ( 1 < il )

    ipnt = ipntp;
    ipntp = ipntp + il;
    il = floor ( il / 2 );
    ndiv = ndiv * 2;
    ihaf = ipntp;

    for iful = ipnt+2 : 2 : ipntp
      ihaf = ihaf + 1;
      rhs(ihaf+1) = rhs(iful+1) - a_cr(3,iful) * rhs(iful) ...
        - a_cr(1,iful+1) * rhs(iful+2);
    end

  end

  rhs(ihaf+1) = rhs(ihaf+1) * a_cr(2,ihaf+1);
  ipnt = ipntp;

  while ( 0 < ipnt )

    ipntp = ipnt;
    ndiv = floor ( ndiv / 2 );
    il = floor ( n / ndiv );
    ipnt = ipnt - il;
    ihaf = ipntp;

    for ifulm = ipnt+1 : 2 : ipntp
      iful = ifulm + 1;
      ihaf = ihaf + 1;
      rhs(iful+1) = rhs(ihaf+1);
      rhs(ifulm+1) = a_cr(2,ifulm+1) * ( rhs(ifulm+1) - a_cr(3,ifulm) * rhs(ifulm) ...
        - a_cr(1,ifulm+1) * rhs(iful+1) );
    end

  end

  x(1:n) = rhs(2:n+1);

  return
end


text:
function [k,c,b]=waitForValidButtonPress

% function [k]=waitForValidButtonPress
% ------------------------------------------------------------------------
% This function is based on the sub-function WFBP found in the |ginput|
% function. It is equivalent to the |waitforbuttonpress| function with the
% exection that that ^C is reserved to end the function as usual.
%
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 03/10/2014
%------------------------------------------------------------------------

% Disabling ^C for edit menu so the only ^C is for interrupting the function.
fig = gcf;
waserr=0;
try

    h=findall(fig,'type','uimenu','accel','C');
    set(h,'accel','');
    keyDown = waitforbuttonpress;
    current_char = double(get(fig,'CurrentCharacter')); % Capturing the character.
    if~isempty(current_char) && (keyDown == 1)           % If the character was generated by the
        if(current_char == 3)                       % current keypress AND is ^C, set 'waserr'to 1
            waserr = 1;                             % so that it errors out.
        end
    end
catch
    waserr = 1;
end
set(h,'accel','C');                                 % Set back the accelerator for edit menu.
drawnow;

if(waserr == 1)
    error('MATLAB:waitForValidButtonPress:Interrupted', 'Interrupted');
end

k = keyDown;
c = get(fig,'CurrentCharacter');
b = abs(get(fig,'CurrentCharacter'));
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function F = coeffs2chebfun2v(X, Y)
%COEFFS2DISKFUNV   Convert componentwise matrices of bivariate Chebyshev
%   coefficients to a chebfun2v.
%
%   F = coeffs2chebfun2v(X, Y) returns a chebfun2v object F that has a
%   bivariate Chebyshev matrices of coefficients X, Y, for each component.
%
% See also CHEBFUN2V/COEFFS2.

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

U = chebfun2(X, 'coeffs');
V = chebfun2(Y, 'coeffs');
F = chebfun2v(U, V);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function mask = compute_conv2_valid_mask(imsz, fsz)

% compute a mask of size IMSZ with 1 entry at valid pixel positions and 0
% invalid oens
%
%   Author:  Stefan Roth, Department of Computer Science, TU Darmstadt
%   Contact: sroth@cs.tu-darmstadt.de
%   $Date:  $
%   $Revision: $

% Copyright 2004-2007, Brown University, Providence, RI. USA
% Copyright 2007-2010 TU Darmstadt, Darmstadt, Germany.
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

mask = ones(imsz);

h1   = floor(fsz(1)/2);
if fsz(1) == 1
    ;
elseif fsz(1) == 2
    mask(end,:) = 0;
elseif mod(fsz(1), 2) == 0
    mask(1:h1-1, :)         = 0;
    mask(end-h1+1:end, :)   = 0;
else
    mask(1:h1, :)    = 0;
    mask(end-h1+1:end, :)   = 0;
end;

h2   = floor(fsz(2)/2);

if fsz(2) == 1
    ;
elseif fsz(2) == 2
    mask(:, end) = 0;
elseif mod(fsz(2), 2) == 0
    mask(:, 1:h2-1)    = 0;
    mask(:, end-h2+1:end)   = 0;
else
    mask(:, 1:h2)    = 0;
    mask(:, end-h2+1:end)   = 0;
end;

author: kristinbranson

repo: JAABA


text:
function Volume=polygon2voxel_double(FacesA,FacesB,FacesC,VerticesX,VerticesY,VerticesZ,VolumeSize,Wrap)
Vertices=[VerticesX(:) VerticesY(:) VerticesZ(:)]-1;

% List with all vertices coordinates of a face
FaceVertices=[Vertices(FacesA,:) Vertices(FacesB,:) Vertices(FacesC,:)];
Volume=false(VolumeSize);
Volume=DrawSplitFaces(FaceVertices,Volume,Wrap);

function Volume=DrawSplitFaces(FaceVertices,Volume,Wrap)
VolumeSize=size(Volume);
% Calculate squared edge distances
dist1=(FaceVertices(:,1)-FaceVertices(:,4)).*(FaceVertices(:,1)-FaceVertices(:,4))+(FaceVertices(:,2)-FaceVertices(:,5)).*(FaceVertices(:,2)-FaceVertices(:,5))+(FaceVertices(:,3)-FaceVertices(:,6)).*(FaceVertices(:,3)-FaceVertices(:,6));
dist2=(FaceVertices(:,7)-FaceVertices(:,4)).*(FaceVertices(:,7)-FaceVertices(:,4))+(FaceVertices(:,8)-FaceVertices(:,5)).*(FaceVertices(:,8)-FaceVertices(:,5))+(FaceVertices(:,9)-FaceVertices(:,6)).*(FaceVertices(:,9)-FaceVertices(:,6));
dist3=(FaceVertices(:,1)-FaceVertices(:,7)).*(FaceVertices(:,1)-FaceVertices(:,7))+(FaceVertices(:,2)-FaceVertices(:,8)).*(FaceVertices(:,2)-FaceVertices(:,8))+(FaceVertices(:,3)-FaceVertices(:,9)).*(FaceVertices(:,3)-FaceVertices(:,9));

% Calculate mFaceVertices(:,1) distance
maxdist=max([dist1(:) dist2(:),dist3(:)],[],2);

% Draw triangle if distance <=0.5 pixel
check=maxdist>0.5;

FVR=FaceVertices(~check,:);
% Select Vertices which must be split
FaceVertices=FaceVertices(check,:);
if(~isempty(FaceVertices))
    dist1=dist1(check);
    dist2=dist2(check);
    dist3=dist3(check);

    DX=(FaceVertices(:,1)+FaceVertices(:,4))/2; DY=(FaceVertices(:,2)+FaceVertices(:,5))/2; DZ=(FaceVertices(:,3)+FaceVertices(:,6))/2;
    FA1=[DX,DY,DZ,FaceVertices(:,4),FaceVertices(:,5),FaceVertices(:,6),FaceVertices(:,7),FaceVertices(:,8),FaceVertices(:,9)];
    FB1=[FaceVertices(:,1),FaceVertices(:,2),FaceVertices(:,3),DX,DY,DZ,FaceVertices(:,7),FaceVertices(:,8),FaceVertices(:,9)];

    DX=(FaceVertices(:,1)+FaceVertices(:,7))/2; DY=(FaceVertices(:,2)+FaceVertices(:,8))/2; DZ=(FaceVertices(:,3)+FaceVertices(:,9))/2;
    FA2=[DX,DY,DZ,FaceVertices(:,4),FaceVertices(:,5),FaceVertices(:,6),FaceVertices(:,7),FaceVertices(:,8),FaceVertices(:,9)];
    FB2=[FaceVertices(:,1),FaceVertices(:,2),FaceVertices(:,3),FaceVertices(:,4),FaceVertices(:,5),FaceVertices(:,6),DX,DY,DZ];

    DX=(FaceVertices(:,7)+FaceVertices(:,4))/2; DY=(FaceVertices(:,8)+FaceVertices(:,5))/2; DZ=(FaceVertices(:,9)+FaceVertices(:,6))/2;
    FA3=[FaceVertices(:,1),FaceVertices(:,2),FaceVertices(:,3),DX,DY,DZ,FaceVertices(:,7),FaceVertices(:,8),FaceVertices(:,9)];
    FB3=[FaceVertices(:,1),FaceVertices(:,2),FaceVertices(:,3),FaceVertices(:,4),FaceVertices(:,5),FaceVertices(:,6),DX,DY,DZ];

    DX=(FaceVertices(:,1)+FaceVertices(:,7))/2; DY=(FaceVertices(:,2)+FaceVertices(:,8))/2; DZ=(FaceVertices(:,3)+FaceVertices(:,9))/2;
    FA4=[DX,DY,DZ,FaceVertices(:,4),FaceVertices(:,5),FaceVertices(:,6),FaceVertices(:,7),FaceVertices(:,8),FaceVertices(:,9)];
    FB4=[FaceVertices(:,1),FaceVertices(:,2),FaceVertices(:,3),FaceVertices(:,4),FaceVertices(:,5),FaceVertices(:,6),DX,DY,DZ];

    dist12=dist1>dist2;
    dist12n=~dist12;
    FA1(dist12n,:)=FA3(dist12n,:);
    FB1(dist12n,:)=FB3(dist12n,:);
    FA2(dist12n,:)=FA4(dist12n,:);
    FB2(dist12n,:)=FB4(dist12n,:);
    dist1(dist12n,:)=dist2(dist12n,:);

    dist13=dist1>dist3;
    dist13n=~dist13;
    FA1(dist13n,:)=FA2(dist13n,:);
    FB1(dist13n,:)=FB2(dist13n,:);

    FaceVertices=[FA1;FB1];

    % Split / Draw Vertices
    Volume=DrawSplitFaces(FaceVertices,Volume,Wrap);
end

% Draw remaining faces
FaceVertices=FVR;

if(Wrap==0)
    % Calculate 1D volume indices
    indexA=mindex3(floor(FaceVertices(:,1)+0.5),floor(FaceVertices(:,2)+0.5), floor(FaceVertices(:,3)+0.5), VolumeSize(1), VolumeSize(2), VolumeSize(3), Wrap);
    indexB=mindex3(floor(FaceVertices(:,4)+0.5),floor(FaceVertices(:,5)+0.5), floor(FaceVertices(:,6)+0.5), VolumeSize(1), VolumeSize(2), VolumeSize(3), Wrap);
    indexC=mindex3(floor(FaceVertices(:,7)+0.5),floor(FaceVertices(:,8)+0.5), floor(FaceVertices(:,9)+0.5), VolumeSize(1), VolumeSize(2), VolumeSize(3), Wrap);

    % Remove outside vertices
    checkA=(FaceVertices(:,1)<0)|(FaceVertices(:,2)<0)|(FaceVertices(:,3)<0)|(FaceVertices(:,1)>(VolumeSize(1)-1))|(FaceVertices(:,2)>(VolumeSize(2)-1))|(FaceVertices(:,3)>(VolumeSize(3)-1));
    checkB=(FaceVertices(:,4)<0)|(FaceVertices(:,5)<0)|(FaceVertices(:,6)<0)|(FaceVertices(:,4)>(VolumeSize(1)-1))|(FaceVertices(:,5)>(VolumeSize(2)-1))|(FaceVertices(:,6)>(VolumeSize(3)-1));
    checkC=(FaceVertices(:,7)<0)|(FaceVertices(:,8)<0)|(FaceVertices(:,9)<0)|(FaceVertices(:,7)>(VolumeSize(1)-1))|(FaceVertices(:,8)>(VolumeSize(2)-1))|(FaceVertices(:,9)>(VolumeSize(3)-1));
    indexA(checkA)=[];
    indexB(checkB)=[];
    indexC(checkC)=[];

    % Draw the vertices
    Volume(indexA)=true;
    Volume(indexB)=true;
    Volume(indexC)=true;
else
     Volume(mindex3(floor(FaceVertices(:,1)+0.5),floor(FaceVertices(:,2)+0.5), floor(FaceVertices(:,3)+0.5), VolumeSize(1), VolumeSize(2), VolumeSize(3), Wrap))=1;
     Volume(mindex3(floor(FaceVertices(:,4)+0.5),floor(FaceVertices(:,5)+0.5), floor(FaceVertices(:,6)+0.5), VolumeSize(1), VolumeSize(2), VolumeSize(3), Wrap))=1;
     Volume(mindex3(floor(FaceVertices(:,7)+0.5),floor(FaceVertices(:,8)+0.5), floor(FaceVertices(:,9)+0.5), VolumeSize(1), VolumeSize(2), VolumeSize(3), Wrap))=1;
end

function index=mindex3(x, y, z, sizx, sizy, sizz, Wrap)
if(Wrap==1)
    % Positive modules
    x=mod(x,sizx);
    y=mod(y,sizy);
    z=mod(z,sizz);
elseif(Wrap>1)
    % Clamp
    x=max(x,0); x=min(x,sizx-1);
    y=max(y,0); y=min(y,sizy-1);
    z=max(z,0); z=min(z,sizz-1);
end
index=z*sizx*sizy+y*sizx+x;
% matlab
index=index+1;

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function y = idualtree3D(w, J, Fsf, sf)

% Inverse 3D Dual-Tree Discrete Wavelet Transform
%
% USAGE:
%   y = idualtree3D(w, J, Fsf, sf)
% INPUT:
%   w - wavelet coefficients
%   J - number of stages
%   Fsf - synthesis filter for the last stage
%   sf - synthesis filters for the preceeding stages
% OUTPUT:
%   y - output arry
% See dualtree3D
%
% WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
% http://taco.poly.edu/WaveletSoftware/

for k = 1:J
    for m = 1:7
        [w{k}{1}{m} w{k}{2}{m} w{k}{3}{m} w{k}{4}{m}] = ...
            pm4inv(w{k}{1}{m}, w{k}{2}{m}, w{k}{3}{m}, w{k}{4}{m});
    end
end

M = [
    1 1 1
    2 2 1
    2 1 2
    1 2 2
];

% initialize output array
y = zeros(2^J * size(w{J}{1}{1}));

for i = 1:4
    f1 = M(i,1);
    f2 = M(i,2);
    f3 = M(i,3);
    yi = w{J+1}{i};
    for k = J:-1:2
        yi = sfb3D(yi, w{k}{i}, sf{f1}, sf{f2}, sf{f3});
    end
    yi = sfb3D(yi, w{1}{i}, Fsf{f1}, Fsf{f2}, Fsf{f3});
    y = y + yi;
end

% normalization
y = y/2;

author: thomaskuestner


text:
function varargout = imRAG(img, varargin)
%IMRAG Region adjacency graph of a label image
%
%   Usage:
%   ADJ = imRAG(LBLMAP);
%   [NODES, ADJ] = imRAG(LBLMAP);
%   [NODES, ADJ, BNDINDS] = imRAG(LBLMAP);
%
%   ADJ = imRAG(LBLMAP);
%   computes region adjacency graph of the 2D or 3D image LBLMAP.
%   The result is a N-by-2 array, containing 2 indices for each couple of
%   neighbor regions. Two regions are considered as neighbor if they are
%   separated by a black (i. e. with color 0) pixel in the horizontal or
%   vertical direction.
%   ADJ has the format [LBL1 LBL2], LBL1 and LBL2 being vertical arrays the
%   same size.
%
%   LBL1 is given in ascending order, LBL2 is given in ascending order for
%   each LBL1. Ex:
%   [1 2]
%   [1 3]
%   [1 4]
%   [2 3]
%   [2 5]
%   [3 4]
%
%   [NODES, ADJ] = imRAG(LBLMAP);
%   Return two arrays: the first one is a [N-by-2] array containing
%   centroids of the N label region, and ADJ is the adjacency list as
%   previously described.
%   For 3D images, the nodes array is [N-by-3].
%
%   [NODES, ADJ, BNDINDS] = imRAG(LBLMAP);
%   (only for 2D arrays) Also returns for each pair of adjacent regions the
%   linear indices of the boundary pixels or voxels (located in between the
%   two regions). BNDINDS is a cell array with as many cells as the number
%   of rows of the ADJ result.
%
%   ... = imRAG(LBLMAP, GAP);
%   Specifies the gap to use between regions. Default is 1, corresponding
%   to regions connected when they are separating by a single 0 pixel or
%   voxel.
%   Using a GAP equal to 0 may be useful when there is no separation
%   between regions. In this case, the result corresponds to the detection
%   of 4-adjacent regions.
%   Using values of GAP larger than 1 allows for detecting adjacent regions
%   when they are separated by more space. However, some adjacencies may be
%   missed in case of thin structures (See algorithm below).
%
%   Example
%     % Simple example on 2D array
%     img = zeros([7 7]);
%     img(2:3, 2:3) = 1; img(5:6, 2:3) = 2;
%     img(2:3, 5:6) = 3; img(5:6, 5:6) = 4;
%     img
%     img =
%          0     0     0     0     0     0     0
%          0     1     1     0     3     3     0
%          0     1     1     0     3     3     0
%          0     0     0     0     0     0     0
%          0     2     2     0     4     4     0
%          0     2     2     0     4     4     0
%          0     0     0     0     0     0     0
%     imRAG(img)
%     ans =
%          1     2
%          1     3
%          2     4
%          3     4
%
%     % Compute the Skeleton of Influence Zone of disjoint regions
%     % (requires image processing toolbox)
%     % read and display an image with several objects
%     img = imread('coins.png');
%     figure(1); clf;
%     imshow(img); hold on;
%     % compute the Skeleton by influence zones using watershed
%     bin = imfill(img>100, 'holes');
%     dist = bwdist(bin);
%     wat = watershed(dist, 4);
%     % compute overlay image for display
%     tmp = uint8(double(img).*(wat>0));
%     ovr = uint8(cat(3, max(img, uint8(255*(wat==0))), tmp, tmp));
%     imshow(ovr);
%     % show the resulting graph
%     [n, e] = imRAG(wat);
%     for i = 1:size(e, 1)
%         plot(n(e(i,:), 1), n(e(i,:), 2), 'linewidth', 4, 'color', 'g');
%     end
%     plot(n(:,1), n(:,2), 'bo', 'markerfacecolor', 'b');
%
%
%   % Create a basic 3D image with labels, and compute RAG
%     germs = [50 50 50;...
%         20 20 20;80 20 20;20 80 20;80 80 20; ...
%         20 20 80;80 20 80;20 80 80;80 80 80];
%     img = zeros([100 100 100]);
%     for i = 1:size(germs, 1)
%         img(germs(i,1), germs(i,2), germs(i,3)) = 1;
%     end
%     wat = watershed(bwdist(img), 6);
%     [n, e] = imRAG(wat);
%     figure; drawGraph(n, e);
%     view(3);
%
%
%   Algorithms
%   The adjacencies are computed by computing the pairs of pixels or voxels
%   with different values and separated by a given distance in each main
%   direction of the image. When values are different, they areassumed to
%   correspond to the labels of two adjacent regions.
%   Adjacencies are detected in the X and Y directions for 3D images, and
%   in the X, Y and Z directions for 3D images. There is no detection of
%   adjacency in diagonal directions.
%
%   An alternative algorithm would be to compute the result of dilation for
%   each region, and detect which labels intersect the dilated region.
%   The number of image processing operations may be very large, ususally
%   resulting in much slower processing time. The following code computes
%   RAG using a 5-by-5 square for detecting adjacent regions:
%     inds = unique(img)
%     inds(inds == 0) = [];
%     adjList = zeros(0, 2);
%     for i = 1:length(inds)
%         mask = imdilate(img == inds(i), ones(5, 5));
%         neighs = unique(img(mask));
%         neighs = neighs(~ismember(neighs, [0 i]));
%         adjList = [adjList; repmat(i, length(neighs), 1) neighs];
%     end
%     adjList = unique(sort(adjList, 2), 'rows');
%   Note that more adjacencies may be detected compared to with the imRAG
%   function.
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2004-02-20,
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

%   History
%   2007-10-12 update doc
%   2007-10-17 add example
%   2010-03-08 replace calls to regionprops by local centroid computation
%   2010-07-29 update doc
%   2012-07-20 remove the use of "diff", using less memory
%   2017-07-31 return output also when nargin == 0, sort labels in 3D

%% Initialisations

% size of image
dim = size(img);

% number of dimensions
nd = length(dim);

% Number of background pixels or voxels between two regions
% gap = 0 -> regions are contiguous
% gap = 1 -> there is a 1-pixel large line or surface between two adjacent
% 	pixels, for example the result of a watershed
gap = 1;
if ~isempty(varargin) && isnumeric(varargin{1})
    gap = varargin{1};
end
shift = gap + 1;

% flag indicating whether edge indices should be computed
computeEdgeInds = nargout > 2;
if computeEdgeInds && gap ~= 1
    error('imRAG:wrongGapValue', ...
        'Edge indices can only be computed for gap equal to 1');
end

if nd == 2
    %% First direction of 2D image

    % identify transitions
    [i1, i2] = find(img(1:end-shift,:) ~= img((shift+1):end, :));

	% get values of consecutive changes
	val1 = img(sub2ind(dim, i1, i2));
	val2 = img(sub2ind(dim, i1+shift, i2));

    % keep only changes not involving background, ordered such that n1 < n2
    inds = val1 ~= 0 & val2 ~= 0 & val1 ~= val2;
    edges = sort([val1(inds) val2(inds)], 2);

    % keep array of positions as linear indices
    if computeEdgeInds
        posD1 = sub2ind(dim, i1(inds)+1, i2(inds));
    end

    %% Second direction of 2D image

    % identify transitions
    [i1, i2] = find(img(:, 1:end-shift) ~= img(:, (shift+1):end));

	% get values of consecutive changes
	val1 = img(sub2ind(dim, i1, i2));
	val2 = img(sub2ind(dim, i1, i2+shift));

    % keep only changes not involving background, ordered such that n1 < n2
    inds = val1 ~= 0 & val2 ~= 0 & val1 ~= val2;
    edges = [edges ; sort([val1(inds) val2(inds)], 2)];

    if computeEdgeInds
        % keep array of positions as linear indices
        posD2 = sub2ind(dim, i1(inds), i2(inds)+1);
        posList = [posD1 ; posD2];
    end

elseif nd == 3
    %% First direction of 3D image

    % identify transitions
    [i1, i2, i3] = ind2sub(dim-[shift 0 0], ...
        find(img(1:end-shift,:,:) ~= img((shift+1):end,:,:)));

	% get values of consecutive changes
	val1 = img(sub2ind(dim, i1, i2, i3));
	val2 = img(sub2ind(dim, i1+shift, i2, i3));

    % keep only changes not involving background
    inds = val1 ~= 0 & val2 ~= 0 & val1 ~= val2;
    edges = unique(sort([val1(inds) val2(inds)], 2), 'rows');

    if computeEdgeInds
        % keep array of positions as linear indices
        posD1 = sub2ind(dim, i1(inds)+1, i2(inds), i3(inds));
    end

    %% Second direction of 3D image

    % identify transitions
    [i1, i2, i3] = ind2sub(dim-[0 shift 0], ...
        find(img(:,1:end-shift,:) ~= img(:,(shift+1):end,:)));

	% get values of consecutive changes
	val1 = img(sub2ind(dim, i1, i2, i3));
	val2 = img(sub2ind(dim, i1, i2+shift, i3));

    % keep only changes not involving background
    inds = val1 ~= 0 & val2 ~= 0 & val1 ~= val2;
    edges = [edges; unique(sort([val1(inds) val2(inds)], 2), 'rows')];

    if computeEdgeInds
        % keep array of positions as linear indices
        posD2 = sub2ind(dim, i1(inds), i2(inds)+1, i3(inds));
    end

    %% Third direction of 3D image

    % identify transitions
    [i1, i2, i3] = ind2sub(dim-[0 0 shift], ...
        find(img(:,:,1:end-shift) ~= img(:,:,(shift+1):end)));

	% get values of consecutive changes
	val1 = img(sub2ind(dim, i1, i2, i3));
    val2 = img(sub2ind(dim, i1, i2, i3+shift));

    % keep only changes not involving background
    inds = val1 ~= 0 & val2 ~= 0 & val1 ~= val2;
    edges = [edges; unique(sort([val1(inds) val2(inds)], 2), 'rows')];

    if computeEdgeInds
        % keep array of positions as linear indices
        posD3 = sub2ind(dim, i1(inds), i2(inds), i3(inds)+1);
        posList = [posD1 ; posD2 ; posD3];
    end
end

% remove double edges, keeping in indsC indices of merged edge for each
% original edge
[edges, indsA, indsC] = unique(edges, 'rows'); %#ok<ASGLU>

if computeEdgeInds
    nEdges = size(edges, 1);
    edgeInds = cell(nEdges, 1);
    for iEdge = 1:nEdges
        inds = indsC == iEdge;
        edgeInds{iEdge} = unique(posList(inds));
    end
end

%% Output processing

if nargout <= 1
    varargout{1} = edges;

else
    % Also compute region centroids
    N = max(img(:));
    points = zeros(N, nd);
    labels = unique(img);
    labels(labels==0) = [];

    if nd == 2
        % compute 2D centroids
        for i = 1:length(labels)
            label = labels(i);
            [iy, ix] = ind2sub(dim, find(img==label));
            points(label, 1) = mean(ix);
            points(label, 2) = mean(iy);
        end
    else
        % compute 3D centroids
        for i = 1:length(labels)
            label = labels(i);
            [iy, ix, iz] = ind2sub(dim, find(img==label));
            points(label, 1) = mean(ix);
            points(label, 2) = mean(iy);
            points(label, 3) = mean(iz);
        end
    end

    % setup output arguments
    varargout{1} = points;
    varargout{2} = edges;

    % eventually returns the position of edges as third output argument
    if nargout > 2
        varargout{3} = edgeInds;
    end
end

author: radishgiant

repo: ThresholdAndSegment


text:
clear; close all; clc;

%cfig = figure('Position', [10,10,1280,1080]);
cfig = figure(1);

% Lidar parameters
lidar = SetLidarParameters();

% Map parameters
borderSize      = 1;            % m
pixelSize       = 0.2;          % m
miniUpdated     = false;        %
miniUpdateDT    = 0.1;          % m
miniUpdateDR    = deg2rad(5);   % rad
% If the robot has moved 0.1 m or rotated 5 degree from last key scan,
% we would add a new key scan and update the map

% Scan matching parameters
fastResolution  = [0.05; 0.05; deg2rad(0.5)]; % [m; m; rad]
bruteResolution = [0.01; 0.01; deg2rad(0.1)]; % not used

% Load lidar data
lidar_data = load('dataset/horizental_lidar.mat');
N = size(lidar_data.timestamps, 1);

% Create an empty map
map.points = [];
map.connections = [];
map.keyscans = [];
pose = [0; 0; 0];
path = pose;

% Here we go!!!!!!!!!!!!!!!!!!!!
for scanIdx = 1 : 1 : N

    disp(['scan ', num2str(scanIdx)]);

    % Get current scan [x1,y1; x2,y2; ...]
    time = lidar_data.timestamps(scanIdx) * 1e-9;
    scan = ReadAScan(lidar_data, scanIdx, lidar, 24);

    % If it's the first scan, initiate
    if scanIdx == 1
        map = Initialize(map, pose, scan);
        miniUpdated = true;
        continue;
    end

    % ===== Matching current scan to local map ============
    % 1. If we executed a mini update in last step, we shall update the
    %    local points map and local grid map (coarse)
    if miniUpdated
        localMap = ExtractLocalMap(map.points, pose, scan, borderSize);
        gridMap1 = OccuGrid(localMap, pixelSize);
        gridMap2 = OccuGrid(localMap, pixelSize/2);
    end

    % 2. Predict current pose using constant velocity motion model
    if scanIdx > 2
        pose_guess = pose + DiffPose(path(:,end-1), pose);
    else
        pose_guess = pose;
    end

    % 3. Fast matching
    if miniUpdated
        [pose, ~] = FastMatch(gridMap1, scan, pose_guess, fastResolution);
    else
        [pose, ~] = FastMatch(gridMap2, scan, pose_guess, fastResolution);
    end

    % 4. Refine the pose using smaller pixels
    % gridMap = OccuGrid(localMap, pixelSize/2);
    [pose, hits] = FastMatch(gridMap2, scan, pose, fastResolution/2);
    %----------------------------------------------------------------------

    % Execute a mini update, if the robot has moved a certain distance
    dp = abs(DiffPose(map.keyscans(end).pose, pose));
    if dp(1)>miniUpdateDT || dp(2)>miniUpdateDT || dp(3)>miniUpdateDR
        miniUpdated = true;
        [map, pose] = AddAKeyScan(map, gridMap2, scan, pose, hits,...
                        pixelSize, bruteResolution, 0.1, deg2rad(3));
    else
        miniUpdated = false;
    end
    path = [path, pose];

    % ===== Loop Closing =========================================
%     if miniUpdated
%         if TryLoopOrNot(map)
%             map.keyscans(end).loopTried = true;
%             map = DetectLoopClosure(map, scan, hits, 4, pi/6, pixelSize);
%         end
%     end

    %----------------------------------------------------------------------

    % Plot
    if mod(scanIdx, 30) == 0
        PlotMap(cfig, map, path, scan, scanIdx);
    end

end

author: meyiao

repo: LaserSLAM

save_path: github-repos/MATLAB/meyiao-LaserSLAM


text:
function J=calcSpherJacob(x,systemType,useHalfRange,lTx,lRx,M)
%%CALCSPHERJACOB Calculate the Jacobian for a monostatic or bistatic
%           spherical measurement with respect to 3D Cartesian position.
%
%INPUTS: x The 3X1 position of the target in Cartesian coordinates in the
%          order [x;y;z].
% systemType An optional parameter specifying the axis from which the
%          angles are measured in radians. Possible values are
%          0 (The default if omitted) Azimuth is measured
%            counterclockwise from the x-axis in the x-y plane. Elevation
%            is measured up from the x-y plane (towards the z-axis). This
%            is consistent with common spherical coordinate systems for
%            specifying longitude (azimuth) and geocentric latitude
%            (elevation).
%          1 Azimuth is measured counterclockwise from the z-axis in the
%            z-x plane. Elevation is measured up from the z-x plane
%            (towards the y-axis). This is consistent with some spherical
%            coordinate systems that use the z axis as the boresight
%            direction of the radar.
%          2 This is the same as 0 except instead of being given
%            elevation, one desires the angle away from the z-axis, which
%            is (pi/2-elevation).
%          3 This is the same as 0 except azimuth is measured clockwise
%            from the y-axis in the x-y plane instead of counterclockwise
%            from the x-axis. This coordinate system often arises when
%            given "bearings" in a local East-North-Up coordinate system,
%            where the bearing directions are measured East of North.
% useHalfRange An optional boolean value specifying whether the bistatic
%          (round-trip) range value has been divided by two. This normally
%          comes up when operating in monostatic mode (the most common
%          type of spherical coordinate system), so that the range
%          reported is a one-way range (or just half a bistatic range).
%          The default if this parameter is not provided is false if lTx
%          is provided and true if it is omitted (monostatic).
%      lTx The 3X1 transmitter position in the global coordinate system
%          with [x;y;z] components. If omitted or an empty matrix is
%          passed, then a vector of zeros is used.
%      lRx The 3X1 receiver position in the global coordinate system
%          with [x;y;z] components. If omitted or an empty matrix is
%          passed, then a vector of zeros is used.
%        M A 3X3 rotation matrices to go from the alignment of the global
%          coordinate system to that at the receiver. If omitted, then it
%          is assumed that the local coordinate system is aligned with the
%          global and M=eye(3) --the identity matrix is used.
%
%OUTPUTS: J A 3X3 Jacobian matrix where the rows are
%          [bistatic range;azimuth;elevation] in that order and the columns
%          take the derivative of the row component with respect to
%          [x,y,z] in that order.
%
%This function just calls rangeGradient and spherAngGradient.
%
%February 2017 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<6||isempty(M))
    M=eye(3,3);
end

if((nargin<4||isempty(lTx))&&(nargin<3||isempty(useHalfRange)))
    useHalfRange=true;
elseif(nargin<3||isempty(useHalfRange))
    useHalfRange=false;
end

if(nargin<5||isempty(lRx))
    lRx=zeros(3,1);
end

if(nargin<4||isempty(lTx))
    lTx=zeros(3,1);
end

if(nargin<2||isempty(systemType))
    systemType=0;
end

J=[rangeGradient(x(1:3,:),useHalfRange,lTx(1:3),lRx(1:3));
   spherAngGradient(x(1:3,:),systemType,lRx(1:3),M)];
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function output = base64decode(varargin)
%
% output = base64decode(input)
%
% Decoding a Base64-encoded byte-stream to recover the original data
% This function depends on JVM in MATLAB or, can optionally use the ZMat
% toolbox (http://github.com/fangq/zmat)
%
% Copyright (c) 2012, Kota Yamaguchi
% URL: https://www.mathworks.com/matlabcentral/fileexchange/39526-byte-encoding-utilities
%
% Modified by: Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%      input: a base64-encoded string
%
% output:
%      output: the decoded binary byte-stream as a uint8 vector
%
% examples:
%      bytes=base64encode('Test JSONLab');
%      orig=char(base64decode(bytes))
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

if(nargin==0)
    error('you must provide at least 1 input');
end
if(exist('zmat','file')==2 || exist('zmat','file')==3)
    output=zmat(varargin{1},0,'base64');
    return;
elseif(isoctavemesh)
    error('You must install the ZMat toolbox (http://github.com/fangq/zmat) to use this function in Octave');
end

error(javachk('jvm'));

if(ischar(varargin{1}))
    varargin{1}=uint8(varargin{1});
end

input=typecast(varargin{1}(:)','uint8');

output = typecast(org.apache.commons.codec.binary.Base64.decodeBase64(input), 'uint8')';

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function output = callpensdpm(interfacedata);

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
ub      = interfacedata.ub;
lb      = interfacedata.lb;

% Bounded variables converted to constraints
if ~isempty(ub)
    [F_struc,K] = addStructureBounds(F_struc,K,ub,lb);
end

% *******************************
% CONVERT FROM INTERNAL FORMAT
% *******************************
penstruct = sedumi2pen(F_struc,K,c,x0);

% **************************
% COPY OPTIONS
% **************************
ops = struct2cell(options.pensdp);ops = [ops{1:end}];
penstruct.ioptions = ops(1:8);
penstruct.foptions = ops(9:end);
penstruct.ioptions(4) = max(0,min(3,options.verbose+1));
if penstruct.ioptions(4)==1
    penstruct.ioptions(4)=0;
end

% FIX
if penstruct.mconstr == 0
    penstruct.msizes = [];
end

if options.savedebug
    save pensdpdebug penstruct
end

%**************************
% CALL PENSDP
%**************************
showprogress('Calling PENSDP',options.showprogress);
solvertime = tic;
[x, fx, u, iresults, fresults, iflag] = pen(penstruct);
solvertime = toc(solvertime);

% Get dual variable (this must be possible to do easier...)
u = u(:);
D_struc = u(1:1:K.l);
if length(K.s)>0
    if K.s(1)>0
        pos = K.l+1;
        for i = 1:length(K.s)
            temp = zeros(K.s(i),K.s(i));
            vecZ = u(pos:pos+0.5*K.s(i)*(K.s(i)+1)-1);
            top = 1;
            for j = 1:K.s(i)
                len = K.s(i)-j+1;
                temp(j:end,j)=vecZ(top:top+len-1);top=top+len;
            end
            temp = (temp+temp');j = find(speye(K.s(i)));temp(j)=temp(j)/2;
            D_struc = [D_struc;temp(:)];
            pos = pos + (K.s(i)+1)*K.s(i)/2;
        end
    end
end

switch iflag
case 0
    problem = 0;
case 1
    problem = 4;
case 2
    problem = 1;
case 3
    problem = 4;
case 4
    problem = 3;
case 5
    problem = 7;
case 6
    problem = 11;
case 7
    problem = 9;
otherwise
    problem = -1;
end
infostr = yalmiperror(problem,'PENSDP/TOMLAB');

if options.savesolveroutput
    solveroutput.f = f;
    solveroutput.x = x;
    solveroutput.u = u;
    solveroutput.iflag = iflag;
    solveroutput.niter = niter;
    solveroutput.feas = feas;
else
    solveroutput = [];
end
if options.savesolverinput
    solverinput.penstruct = penstruct;
else
    solverinput = [];
end

% Standard interface
output = createOutputStructure(x(:),D_struc,[],problem,infostr,solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function res = isCounterClockwise(p1, p2, p3, varargin)
%ISCOUNTERCLOCKWISE Compute relative orientation of 3 points
%
%   CCW = isCounterClockwise(P1, P2, P3);
%   Computes the orientation of the 3 points. The returns is:
%   +1 if the path P1->P2->P3 turns Counter-Clockwise (i.e., the point P3
%       is located "on the left" of the line P1-P2)
%   -1 if the path turns Clockwise (i.e., the point P3 lies "on the right"
%       of the line P1-P2)
%   0  if the point P3 is located on the line segment [P1 P2].
%
%   This function can be used in more complicated algorithms: detection of
%   line segment intersections, convex hulls, point in triangle...
%
%   CCW = isCounterClockwise(P1, P2, P3, EPS);
%   Specifies the threshold used for detecting colinearity of the 3 points.
%   Default value is 1e-12 (absolute).
%
%   Example
%   isCounterClockwise([0 0], [10 0], [10 10])
%   ans =
%       1
%   isCounterClockwise([0 0], [0 10], [10 10])
%   ans =
%       -1
%   isCounterClockwise([0 0], [10 0], [5 0])
%   ans =
%       0
%
%   See also
%   points2d, isPointOnLine, isPointInTriangle
%
%   References
%   Algorithm adapated from Sedgewick's book.
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-04-09
% Copyright 2011 INRA - Cepia Software Platform.

%   HISTORY
%   2011-05-16 change variable names, add support for point arrays

% get threshold value
eps = 1e-12;
if ~isempty(varargin)
    eps = varargin{1};
end

% ensure all data have same size
np = max([size(p1, 1) size(p2, 1) size(p3,1)]);
if np > 1
    if size(p1,1) == 1
        p1 = repmat(p1, np, 1);
    end
    if size(p2,1) == 1
        p2 = repmat(p2, np, 1);
    end
    if size(p3,1) == 1
        p3 = repmat(p3, np, 1);
    end
end

% init with 0
res = zeros(np, 1);

% extract vector coordinates
x0  = p1(:, 1);
y0  = p1(:, 2);
dx1 = p2(:, 1) - x0;
dy1 = p2(:, 2) - y0;
dx2 = p3(:, 1) - x0;
dy2 = p3(:, 2) - y0;

% check non colinear cases
res(dx1 .* dy2 > dy1 .* dx2) =  1;
res(dx1 .* dy2 < dy1 .* dx2) = -1;

% case of colinear points
ind = abs(dx1 .* dy2 - dy1 .* dx2) < eps;
res(ind( (dx1(ind) .* dx2(ind) < 0) | (dy1(ind) .* dy2(ind) < 0) )) = -1;
res(ind(  hypot(dx1(ind), dy1(ind)) <  hypot(dx2(ind), dy2(ind)) )) =  1;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function test1 (nmat)
%TEST1 test for BTF
% Requires CSparse and UFget
% Example:
%   test1
% See also btf, maxtrans, strongcomp, dmperm, UFget,
%   test1, test2, test3, test4, test5.

% Copyright 2007, Timothy A. Davis, University of Florida

index = UFget ;
% f = find (index.sprank < min (index.nrows, index.ncols)) ;
f = 1:length (index.nrows) ;

% too much time:
skip = [1514 1297 1876 1301] ;

f = setdiff (f, skip) ;

[ignore i] = sort (index.nnz (f)) ;
f = f (i) ;

if (nargin < 1)
    nmat = 1000 ;
end
nmat = min (nmat, length (f)) ;
f = f (1:nmat) ;

T0 = zeros (nmat,1) ;
T1 = zeros (nmat,1) ;
Anz = zeros (nmat,1) ;
figure (1) ;
MN = zeros (nmat, 2) ;
Nzdiag = zeros (nmat,1) ;

% warmup
p = maxtrans (sparse (1)) ;         %#ok
p = cs_dmperm (sparse (1)) ;        %#ok
a = cs_transpose (sparse (1)) ;     %#ok

h = waitbar (0, 'BTF test 1 of 6') ;

try

    for k = 1:nmat

        Prob = UFget (f (k), index) ;
        A = Prob.A ;
        clear Prob
        t = 0 ;

        waitbar (k/nmat, h) ;

        r = full (sum (spones (A), 2)) ;
        c = full (sum (spones (A))) ;
        m2 = length (find (r > 0)) ;
        n2 = length (find (c > 0)) ;

        if (m2 < n2)
            tic
            A = cs_transpose (A) ;
            t = toc ;
        end

        Nzdiag (k) = nnz (diag (A)) ;

        [m n] = size (A) ;
        Anz (k) = nnz (A) ;
        MN (k,:) = [m n] ;

        tic
        q = maxtrans (A) ;
        t0 = toc ;
        s0 = sum (q > 0) ;
        T0 (k) = max (1e-9, t0) ;

        tic
        p = cs_dmperm (A) ;
        t1 = toc ;
        s1 = sum (p > 0) ;
        T1 (k) = max (1e-9, t1) ;

        fprintf (...
            '%4d maxtrans %10.6f %10.6f  cs_dmperm %10.6f m/n %8.2f', ...
            f(k), t, t0, t1, m/n) ;
        if (t1 ~= 0)
            fprintf (' rel: %8.4f', t0 / t1) ;
        end
        fprintf ('\n') ;
        if (s0 ~= s1)
            error ('!') ;
        end

        if (s0 == n & m == n)                                               %#ok
            B = A (:, q) ;
            subplot (2,2,1) ;
            cspy (B) ;
            if (nnz (diag (B)) ~= n)
                error ('?')
            end
            clear B
        else
            cspy (0) ;
        end

        maxnz = nnz (A) ;

        zfree  = find (MN (1:k,1) == MN (1:k,2) & Nzdiag (1:k) == MN(1:k,1)) ;
        square = find (MN (1:k,1) == MN (1:k,2) & Nzdiag (1:k) ~= MN(1:k,1)) ;
        tall   = find (MN (1:k,1) >  MN (1:k,2)) ;
        squat  = find (MN (1:k,1) <  MN (1:k,2)) ;

        subplot (2,2,2) ;
        loglog (Anz (square), T0 (square) ./ T1 (square), ...
            'o', [1 maxnz], [1 1], 'r-') ;
        title ('square') ;
        subplot (2,2,3) ;
        loglog (Anz (tall), T0 (tall) ./ T1 (tall), ...
            'o', [1 maxnz], [1 1], 'r-') ;
        title ('tall') ;
        subplot (2,2,4) ;
        title ('square, intially zero-free') ;
        loglog (Anz (zfree), T0 (zfree) ./ T1 (zfree), ...
            'o', [1 maxnz], [1 1], 'r-') ;
        title ('square, zero-free diag') ;

        drawnow

    end

catch
    % out-of-memory is OK, other errors are not
    disp (lasterr) ;
    if (isempty (strfind (lasterr, 'Out of memory')))
        error (lasterr) ;                                                   %#ok
    else
        fprintf ('test terminated early, but otherwise OK\n') ;
    end
end

close (h) ;


text:
function jac = p17_jac ( neqn, t, y )

%*****************************************************************************80
%
%% P17_JAC evaluates the jacobian for problem p17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the jacobian.
%
%    Output, real JAC(NEQN,NEQN), the jacobian matrix.
%
  jac = zeros ( neqn, neqn );

  d = ( sqrt ( ( y(1).^2 + y(2).^2 ) ) ).^5;

  jac(1,3) = 1.0;
  jac(2,4) = 1.0;
  jac(3,1) = ( 2.0 * y(1).^2 - y(2).^2 ) / d;
  jac(3,2) = 3.0 * y(1) * y(2) / d;
  jac(4,1) = 3.0 * y(1) * y(2) / d;
  jac(4,2) = ( - y(1).^2 + 2.0 * y(2).^2 ) / d;

  return
end


text:
function [yhat prdt_per_tree] = classRFPredict( x, rf )
    [nsample ndim] = size(x);
    ntree = rf.ntree;
    prdt_per_tree = zeros(nsample, ntree);

    ntree = rf.ntree;

    for ix = 1 : nsample
        for n = 1 : ntree
            k = 1;
            while rf.nodestatus{n}(k) ~= -1
                if x(ix, rf.bestvar{n}(k)) <= rf.bestsplit{n}(k)
                    k = rf.treemap{n}(k, 1);
                else
                    k = rf.treemap{n}(k, 2);
                end
            end
            prdt_per_tree(ix, n) = rf.nodeclass{n}(k);
            % fprintf( 'ix: %d, n: %d\n', ix, n );
        end
    end

    yhat = zeros(nsample, 1);
    for ix = 1 : nsample
        % for n = 1 : ntree
            % when there are equal votes for the classes, it favors the
            % first class
            yhat(ix) = mode( prdt_per_tree(ix, :) );
        % end
    end

    for ix = 1 : length(rf.new_labels)
        yhat( yhat == rf.new_labels(ix) ) = rf.orig_labels(ix);
    end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm