text:
function SleepState = SleepScoreMaster(basePath,varargin)
%SleepScoreMaster(basePath,<options>)
%This is the master function for sleep state scoring.
%
%It's strongly recommended that you
%   1) indicate bad (noisy) channels using bz_getSessionInfo(basePath,'editGUI',true)
%      before running SleepScoreMaster.
%   2) use the 'ignoretime' input to exclude time windows with opto
%       stimulation or behavior with electrical noise
%   3) check the scoring quality using TheStateEditor after running SleepScoreMaster.
%      Use the 'A' key in TheStateEditory to further refine thresholds as
%      needed, and implement sticky thresholds.
%
%INPUT
%   basePath        folder containing .xml and .lfp files.
%                   basePath and files should be of the form:
%                   'whateverfolder/recordingName/recordingName'
%   (optional)      If no inputs included, select folder(s) containing .lfp
%                   and .xml file in prompt.
%   (optional)      if no .lfp in basePath, option to select multiple
%                   lfp-containing subfolders
%
%   OPTIONS
%   'savedir'       Default: basePath
%   'overwrite'     Overwrite all processing steps (Default: false)
%   'ignoreManual'  Default: false. Overwrite manual scoring from TheStateEditor
%   'savebool'      Default: true. Save anything.
%   'scoretime'     Window of time to score. Default: [0 Inf]
%                   NOTE: must be continous interval
%   'ignoretime'    Time intervals winthin scoretime to ignore
%                   (for example, opto stimulation or behavior with artifacts)
%   'winparms'      [FFT window , smooth window] (Default: [2 15])
%                   (Note: updated from [10 10] based on bimodaility optimization, 6/17/19)
%   'SWWeightsName' Name of file in path (in Dependencies folder)
%                   containing the weights for the various frequencies to
%                   be used for SWS detection.
%                   Default is to use Power Spectrum Slope ('PSS'),
%                   If this doesn't work, try 'SWweights.mat'
%                   or 'SWweightsHPC.mat' for HPC recording.
%   'Notch60Hz'     Boolean 0 or 1.  Value of 1 will notch out the 57.5-62.5 Hz
%                   band, default is 0, no notch.  This can be necessary if
%                   electrical noise.
%   'NotchUnder3Hz' Boolean 0 or 1.  Value of 1 will notch out the 0-3 Hz
%                   band, default is 0, no notch.  This can be necessary
%                   due to poor grounding and low freq movement transients
%   'NotchHVS'      Boolean 0 or 1.  Value of 1 will notch the 4-10 and
%                   12-18 Hz bands for SW detection, default is 0, no
%                   notch.  This can be useful in
%                   recordings with prominent high voltage spindles which
%                   have prominent ~16hz harmonics
%   'NotchTheta'    Boolean 0 or 1.  Value of 1 will notch the 4-10 Hz
%                   band for SW detection, default is 0, no notch.  This
%                   can be useful to
%                   transform the cortical spectrum to approximately
%                   hippocampal, may also be necessary with High Voltage
%                   Spindles
%   'stickytrigger' Implements a "sticky" trigger for SW/EMG threshold
%                   crossings: metrics must reach halfway between threshold
%                   and opposite peak to count as crossing (reduces
%                   flickering) (default:true)
%   'SWChannels'    A vector list of channels that may be chosen for SW
%                   signal
%   'ThetaChannels' A vector list of channels that may be chosen for Theta
%                   signal
%   'rejectChannels' A vector of channels to exclude from the analysis
%   'saveLFP'       (default:true) to save SleepScoreLFP.lfp.mat file
%   'noPrompts'     (default:false) an option to not prompt user of things
%
%OUTPUT
%   !THIS IS OUT OF DATE - UPDATE!
%   StateIntervals  structure containing start/end times (seconds) of
%                   NREM, REM, WAKE states and episodes. states is the
%                   "raw" state scoring. episodes are joined episodes of
%                   extended (40s) time in a given states, allowing for
%                   brief interruptions. also contains NREM packets,
%                   unitary epochs of NREM as described in Watson et al 2016.
%                   saved in a .mat file:
%                   recordingname_SleepScore.mat
%
%
% DLevenstein and BWatson 2015/16
%% Recording Selection
%if recname is 'select' or something
%use uigetfile to pick and get list of filenames
%if recname is 'all', get all recordings in a folder and
%then run SleepScoreMaster on each of the filenames'
%if no input arguements... select uigetfile

%Select from no input
if ~exist('basePath','var')
    basePath = uigetdir(cd,...
        'Which recording(s) would you like to state score?');
    if isequal(basePath,0);return;end
end

%Separate datasetfolder and recordingname
[datasetfolder,recordingname,extension] = fileparts(basePath);
recordingname = [recordingname,extension]; % fileparts parses '.' into extension

%% If there is no .lfp in basePath, choose (multiple?) folders within basePath.
%Select from dataset folder - need to check if .xml/lfp exist
if ~exist(fullfile(datasetfolder,recordingname,[recordingname,'.lfp']),'file') && ...
    ~exist(fullfile(datasetfolder,recordingname,[recordingname,'.eeg']),'file')
    display(['no .lfp file in basePath, pick a selection of session folders',...
             'containing .lfp files'])
         %Find all basePaths within the topPath
        [basePaths,recordingname] = bz_FindBasePaths(basePath,'select',true);
end

%If multiple recordings, loop calling SleepScoreMaster with each
numrecs = length(recordingname);
if numrecs > 1 & iscell(recordingname)
    display(['Multiple Recordings (',num2str(numrecs),')'])
    for rr = 1:numrecs
        multibasepath = basePaths{rr};
        SleepScoreMaster(multibasepath,varargin{:})
        close all
    end
    return
elseif numrecs == 1 & iscell(recordingname)
        recordingname = recordingname{1};
end

display(['Scoring Recording: ',recordingname]);

%% inputParse for Optional Inputs and Defaults
p = inputParser;

addParameter(p,'overwrite',false)
addParameter(p,'savebool',true,@islogical)
addParameter(p,'savedir',datasetfolder)
addParameter(p,'scoretime',[0 Inf])
addParameter(p,'ignoretime',[])
addParameter(p,'SWWeightsName','PSS')
addParameter(p,'Notch60Hz',0)
addParameter(p,'NotchUnder3Hz',0)
addParameter(p,'NotchHVS',0)
addParameter(p,'NotchTheta',0)
addParameter(p,'SWChannels',0)
addParameter(p,'ThetaChannels',0)
addParameter(p,'rejectChannels',[]);
addParameter(p,'noPrompts',true);
addParameter(p,'stickytrigger',false);
addParameter(p,'saveLFP',true);
addParameter(p,'winparms',[2 15]);
addParameter(p,'ignoreManual',false)

parse(p,varargin{:})
%Clean up this junk...
overwrite = p.Results.overwrite;
savedir = p.Results.savedir;
savebool = p.Results.savebool;
scoretime = p.Results.scoretime;
ignoretime = p.Results.ignoretime;
SWWeightsName = p.Results.SWWeightsName;
Notch60Hz = p.Results.Notch60Hz;
NotchUnder3Hz = p.Results.NotchUnder3Hz;
NotchHVS = p.Results.NotchHVS;
NotchTheta = p.Results.NotchTheta;
SWChannels = p.Results.SWChannels;
ThetaChannels = p.Results.ThetaChannels;
rejectChannels = p.Results.rejectChannels;
noPrompts = p.Results.noPrompts;
stickytrigger = p.Results.stickytrigger;
saveLFP = p.Results.saveLFP;
winparms = p.Results.winparms;
ignoreManual = p.Results.ignoreManual;

%% Database File Management
savefolder = fullfile(savedir,recordingname);
if ~exist(savefolder,'dir')
    mkdir(savefolder)
end

%Filenames of metadata and SleepState.states.mat file to save
sessionmetadatapath = fullfile(savefolder,[recordingname,'.SessionMetadata.mat']);
%Buzcode outputs
bz_sleepstatepath = fullfile(savefolder,[recordingname,'.SleepState.states.mat']);

%% Check for existing Manual Scoring
if exist(bz_sleepstatepath,'file') && ~overwrite && ~ignoreManual
   SleepState_old = load(bz_sleepstatepath);
   if isfield(SleepState_old.SleepState.detectorinfo,'LastManualUpdate')
       display(['Manual scoring detected... will update the SleepScoreMetrics and AutoScoreInts, ',...
        'but keep previous (manual) scoring.'])
       display(['To overwrite manual scoring use ''ignoreManual'',true ',...
           'or ''overwrite'',true to overwrite all metrics.'])
       ManScore.ints = SleepState_old.SleepState.ints;
       ManScore.idx = SleepState_old.SleepState.idx;
       ManScore.LastManualUpdate = SleepState_old.SleepState.detectorinfo.LastManualUpdate;
   end
end

%% Get channels not to use
sessionInfo = bz_getSessionInfo(basePath,'noPrompts',noPrompts);
% check that SW/Theta channels exist in rec..
if length(SWChannels) > 1
    if sum(ismember(SWChannels,sessionInfo.channels)) ~= length(SWChannels)
        error('some of the SW input channels dont exist in this recording...?')
    end
end
if length(ThetaChannels) > 1
    if sum(ismember(ThetaChannels,sessionInfo.channels)) ~= length(ThetaChannels)
        error('some of the theta input channels dont exist in this recording...?')
    end
end

%Is this still needed?/Depreciated?
% if exist(sessionmetadatapath,'file')%bad channels is an ascii/text file where all lines below the last blank line are assumed to each have a single entry of a number of a bad channel (base 0)
%     load(sessionmetadatapath)
%     rejectChannels = [rejectChannels SessionMetadata.ExtracellEphys.BadChannels];
% elseif isfield(sessionInfo,'badchannels')
if isfield(sessionInfo,'badchannels')
    rejectChannels = [rejectChannels sessionInfo.badchannels]; %get badchannels from the .xml

end

if isempty(rejectChannels)
    display('No rejected channels - it''s recommended you identify noisy channels to ignore')
end

%% CALCULATE EMG FROM HIGH-FREQUENCY COHERENCE
% Load/Calculate EMG based on cross-shank correlations
% (high frequency correlation signal = high EMG).
% Schomburg E.W. Neuron 84, 470?485. 2014)
EMGFromLFP = bz_EMGFromLFP(basePath,'overwrite',overwrite,...
                                     'rejectChannels',rejectChannels,'noPrompts',noPrompts,...
                                     'saveMat',savebool);

%% DETERMINE BEST SLOW WAVE AND THETA CHANNELS
%Determine the best channels for Slow Wave and Theta separation.
%Described in Watson et al 2016, with modifications
SleepScoreLFP = PickSWTHChannel(basePath,...
                            scoretime,SWWeightsName,...
                            Notch60Hz,NotchUnder3Hz,NotchHVS,NotchTheta,...
                            SWChannels,ThetaChannels,rejectChannels,...
                            overwrite,'ignoretime',ignoretime,...
                            'noPrompts',noPrompts,'saveFiles',saveLFP&savebool,...
                            'window',winparms(1),'smoothfact',winparms(2),'IRASA',true);

%% CLUSTER STATES BASED ON SLOW WAVE, THETA, EMG

%Calculate the scoring metrics: broadbandLFP, theta, EMG
display('Quantifying metrics for state scoring')
[SleepScoreMetrics,StatePlotMaterials] = ClusterStates_GetMetrics(...
                                           basePath,SleepScoreLFP,EMGFromLFP,overwrite,...
                                           'onSticky',stickytrigger,'ignoretime',ignoretime,...
                                           'window',winparms(1),'smoothfact',winparms(2),...
                                           'IRASA',true,'ThIRASA',true);

%Use the calculated scoring metrics to divide time into states
display('Clustering States Based on EMG, SW, and TH LFP channels')
[ints,idx,MinTimeWindowParms] = ClusterStates_DetermineStates(SleepScoreMetrics);

%% RECORD PARAMETERS from scoring
detectionparms.userinputs = p.Results;
detectionparms.MinTimeWindowParms = MinTimeWindowParms;
detectionparms.SleepScoreMetrics = SleepScoreMetrics;

% note and keep special version of original hists and threshs

SleepState.ints = ints;
SleepState.idx = idx;
SleepState.detectorinfo.detectorname = 'SleepScoreMaster';
SleepState.detectorinfo.detectionparms = detectionparms;
SleepState.detectorinfo.detectionparms.histsandthreshs_orig = detectionparms.SleepScoreMetrics.histsandthreshs;
SleepState.detectorinfo.detectiondate = datestr(now,'yyyy-mm-dd');
SleepState.detectorinfo.StatePlotMaterials = StatePlotMaterials;

%Put old manual scoring back in
if exist('ManScore','var')
    SleepState.AutoScoreInts = SleepState.ints;
    SleepState.ints = ManScore.ints;
    SleepState.idx = ManScore.idx;
    SleepState.detectorinfo.LastManualUpdate = ManScore.LastManualUpdate;
end

%Saving SleepStates
save(bz_sleepstatepath,'SleepState');

%% MAKE THE STATE SCORE OUTPUT FIGURE
%ClusterStates_MakeFigure(stateintervals,stateIDX,figloc,SleepScoreMetrics,StatePlotMaterials);
try
    ClusterStates_MakeFigure(SleepState,basePath,noPrompts);
    disp('Figures Saved to StateScoreFigures')
catch
    disp('Figure making error')
end

%% JOIN STATES INTO EPISODES

% Extract states, Episodes, properly organize params etc, prep for final saving
display('Calculating/Saving Episodes')
StatesToEpisodes(SleepState,basePath);

display(['Sleep Score ',recordingname,': Complete!']);

%% PROMPT USER TO MANUALLY CHECK DETECTION WITH THESTATEEDITOR
if ~noPrompts
    str = input('Would you like to check detection with TheStateEditor? [Y/N] ','s');
    switch str
        case {'Y','y',''}
            TheStateEditor([basePath,filesep,recordingname])
        case {'N','n'}
        otherwise
            display('Unknown input..... you''ll have to load TheStateEditor on your own')
    end
end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function roi = dtiReadRoi(fileName, xformParams)
%
%  roi = dtiReadRoi([fileName], [xformParams])
%
% HISTORY:
%  2005.02.24 RFD: pulled code from other places to modularize.
%  2008.04.01 RFD: now forces ROI coords to 'double' type.

if(~exist('fileName','var') || isempty(fileName))
    [f, p] = uigetfile({'*.mat';'*.*'}, 'Load ROI file...');
    if(isnumeric(f)), disp('Read ROI canceled.'); roi=[]; return; end
    fileName = fullfile(p,f);
end
if(~exist('xformParams','var')) xformParams = []; end
load(fileName);

% We now save version and coordinate space fields. But for backwards compatability....
if(~exist('versionNum')) versionNum = 0.1; end
if(~exist('coordinateSpace')) coordinateSpace = 'acpc'; end
if ~exist('roi','var'), error('No ROI variable found.'); end

% New field in ROI struct for the DTI query ID. We alway reset it to -1
% when loading to indicate no associated dtiQuery ROI.
roi.query_id = -1;

roi.coords = double(roi.coords(~any(isnan(roi.coords')),:));
roi.visible = 1;
if(~strcmp(coordinateSpace,'acpc'))
    % try to find a matching coordinate space xform
    if(~isfield(xformParams,'name'))
        if(strcmp(coordinateSpace,'MNI'))
            csInd = 1;
        else
            error('t1NormParams has an old format (no coord space name).');
        end
    else
        csInd = strmatch(coordinateSpace, {xformParams.name}, 'exact');
        if(isempty(csInd))
            warning(['No coord spaces matching ROI coord space "' coordinateSpace '".']);
        end
    end
    if(~isempty(xformParams) && ~isempty(xformParams(csInd)))
        if(~isfield(xformParams(csInd),'sn') || isempty(xformParams(csInd).sn))
            disp(['Subject data appear to be in the same space as the ROI (' coordinateSpace ')- not warping ROI...']);
        else
            disp(['Warping ROI coords from ' coordinateSpace ' space to this subject''s space...']);
            % *** TO DO: we assume that the xform is an SPM-style sn struct. We
            % should figure out how to support other xforms.
            roi.coords = mrAnatXformCoords(xformParams(csInd).sn, roi.coords);
        end
    else
        disp(['NOTE: ROI coordinate space is ' coordinateSpace ', but there are no spatial norm params- ROI loaded WITHOUT coordinate transform.']);
    end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [I, T_ini,T_ref] = lime(L,para,post)
% A wrapper of LIME
% LIME: A Method for Low-light IMage Enhancement
%
% INPUTS
% para: parameter
% post: true if do post-processing (denoise)
% OUTPUTS
% I - result image
% -------------------------------------------------------------------------
% @article{guo2016lime,
%   title={LIME: A Method for Low-light IMage Enhancement},
%   author={Guo, Xiaojie},
%   journal={arXiv preprint arXiv:1605.05034},
%   year={2016}
% }
% -------------------------------------------------------------------------
%
% Zhenqiang Ying
% 2016-11-02

if nargin == 0
   I = imload(uigetimagefile);
   [J, T_ini,T_ref] = lime(I);
   T = repmat(T_ref, [1 1 3]);
   ezFig I J I./(T.^0.6) I./(T.^0.8) T_ini T_ref
   clear I;
   return;
end

Require LIME

if nargin < 2
    para.lambda = 0.15;
    para.sigma = 2;
    para.gamma = 0.8;
    if nargin < 3
       post = false;
    end
end

if ~isfloat(L), L = im2double(L); end

tic
[I, T_ini,T_ref] = LIME(L,para);
toc

% figure(1);imshow(L);title('Input');
% figure(2);imshow(I);title('LIME');

% figure(3);imshow(T_ini,[]);colormap hot;title('Initial T');
% figure(4);imshow(T_ref,[]);colormap hot;title('Refined T');
%% Post Processing
if post
    %addpath(genpath('+existingMethods\+LIME\BM3D'));
    %import filter.BM3D.*

    YUV = rgb2ycbcr(I);
    Y = YUV(:,:,1);

    Y_d = BM3D(Y);

    I_d = ycbcr2rgb(cat(3,Y_d,YUV(:,:,2:3)));
    I_f = (I).*repmat(T_ref,[1,1,3])+I_d.*repmat(1-T_ref,[1,1,3]);

%     figure(5);imshow(I_d);title('Denoised ');
%     figure(6);imshow(I_f);title('Recomposed');
    I = I_f;
end

end

%% DEMO
% clc;close all;clear all;addpath(genpath('./'));
% %%
% filename = '1.bmp';
% L = (im2double(imread(filename)));
%
% post = true;
%
% para.lambda = 0.15;
% para.sigma = 2;
% para.gamma = 0.8;
% tic
% [I, T_ini,T_ref] = LIME(L,para);
% toc
%
% figure(1);imshow(L);title('Input');
% figure(2);imshow(I);title('LIME');
% figure(3);imshow(T_ini,[]);colormap hot;title('Initial T');
% figure(4);imshow(T_ref,[]);colormap hot;title('Refined T');
% %% Post Processing
% if post
% YUV = rgb2ycbcr(I);
% Y = YUV(:,:,1);
%
% sigma_BM3D = 50;
% [~, Y_d] = BM3D(Y,Y,sigma_BM3D,'lc',0);
%
% I_d = ycbcr2rgb(cat(3,Y_d,YUV(:,:,2:3)));
% I_f = (I).*repmat(T_ref,[1,1,3])+I_d.*repmat(1-T_ref,[1,1,3]);
%
% figure(5);imshow(I_d);title('Denoised ');
% figure(6);imshow(I_f);title('Recomposed');
% end

author: baidut

repo: BIMEF

save_path: github-repos/MATLAB/baidut-BIMEF


text:
% This function implements the relaxed sparsification descibed in Section 3.4

function softSegments = sparsifySegments(softSegments, Laplacian, imageGrad)

    sigmaS = 1; % sparsity
    sigmaF = 1; % fidelity
    delta = 100; % constraint
    [h, w, compCnt] = size(softSegments);
    N = h * w * compCnt;

    if ~exist('imageGrad', 'var') || isempty(imageGrad)
        % If image gradient is not provided, set the param to the default 0.9
        spPow = 0.90;
    else
        % Compute the per-pixel sparsity parameter from the gradient
        imageGrad(imageGrad > 0.1) = 0.1;
        imageGrad = imageGrad + 0.9;
        spPow = repmat(imageGrad(:), [compCnt, 1]);
    end

    % Iter count for pcg and main optimization
    itersBetweenUpdate = 100;
    highLevelIters = 20;

    % Get rid of very low/high alpha values and normalize
    softSegments(softSegments < 0.1) = 0;
    softSegments(softSegments > 0.9) = 1;
    softSegments = softSegments ./ repmat(sum(softSegments, 3), [1 1 size(softSegments, 3)]);

    % Construct the linear system
    lap = Laplacian;
    for i = 2 : compCnt
        Laplacian = blkdiag(Laplacian, lap);
    end

    % The alpha constraint
    C = repmat(speye(h*w), [1 compCnt]);
    C = C' * C;
    Laplacian = Laplacian + delta * C;

    % The sparsification optimization
    softSegments = softSegments(:);
    compInit = softSegments; % needed for fidelity energy
    for iter = 1 : highLevelIters
        if rem(iter, 5) == 0
            disp(['               Iteration ' int2str(iter) ' of ' int2str(highLevelIters)]);
        end
        [u, v] = getUandV(softSegments, spPow); % The sparsity energy
        A = Laplacian + sigmaS * (spdiags(u, 0, N, N) + spdiags(v, 0, N, N)) + sigmaF * speye(N);
        b = sigmaS * v + sigmaF * compInit + delta;
        [softSegments, ~] = pcg(A, b, [], itersBetweenUpdate, [], [], softSegments);
    end

    % One final iter for good times (everything but sparsity)
    A = Laplacian + sigmaF * speye(N);
    b = sigmaF * softSegments + delta;
    softSegments = pcg(A, b, [], 10 * itersBetweenUpdate, [], [], softSegments);

    % Ta-dah
    softSegments = reshape(softSegments, [h w compCnt]);
end

function [u, v] = getUandV(comp, spPow)
    % Sparsity terms in the energy
    eps = 1e-2;
    u = max(abs(comp(:)), eps) .^ (spPow - 2);
    v = max(abs(1 - comp(:)), eps) .^ (spPow - 2);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open 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.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

function [robot_i, frame_i] = gtsamFrameIdToIndices(id)
% robot_i and frame_i are 1-indexed!
assert(isa(id, 'int64'));
frame_i = mod(id, 2^56);
robot_i = (id - frame_i + 1) / 2^56 - 'a' + 1;
end

author: uzh-rpg


text:
function ROI = motionCompSaveROI(view,ROI,ROIname,resample)
%
%    gb 05/14/05
%
%    ROI = motionCompSaveROI(view,ROI,ROIname,resample)
%
% Stores the ROI into the directory HOMEDIR/Inplane/ROI with the name
% ROIname.
%
% Inputs:
%   - view: current Inplane view
%   - ROIname: string name of the ROI.
%   - resample: The ROI is stored in the anatomy format which is twice
%   bigger than the functional data format. Set resample to 1 if you are
%   saving a ROI having the size of the functional data.

% Initializes arguments and variables
global dataTYPES
curDataType = viewGet(view,'currentDataType');
curScan = viewGet(view,'currentScan');

if ieNotDefined('ROIname')
    ROIname = 'ROInoname';
end

if ieNotDefined('resample')
    resample = 1;
end

scan = viewGet(view,'curScan');

% Resamples the ROI if needed
if resample
    size1 = sliceDims(view,scan);
    size2 = size(view.anat);
    T = maketform('affine',[size2(2)/size1(2) 0 0; 0 size2(1)/size1(1) 0; 0 0 1]);

    ROI = imtransform(ROI,T,'nearest');
    ROI = [zeros(1,size(ROI,2) + 1,size(ROI,3));zeros(size(ROI,1),1,size(ROI,3)),ROI];

    nVoxels = size2(1:2);
else
    nVoxels = sliceDims(view,scan);
end

% Transforms it to the storage format which is a set of 3-coordinate
% vectors belonging to the ROI.
roiIndex = find(ROI)' - 1;
ROIz = floor(roiIndex/prod(nVoxels));
roiIndex = rem(roiIndex,prod(nVoxels));
ROIy = floor(roiIndex/nVoxels(1));
ROIx = rem(roiIndex,nVoxels(1));

ROIx = ROIx + 1;
ROIy = ROIy + 1;
ROIz = ROIz + 1;
ROI = [ROIx;ROIy;ROIz];

% Saves the ROI
pathStr = roiDir(view);
pathStr = fullfile(pathStr,[ROIname '.mat']);
ROI = struct('color','b','coords',ROI,'name',ROIname,'viewType','Inplane');
save(pathStr,'ROI');


text:
function inspect_bug1230

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

% this is not really pertaining to a bug, but to a feature request.

load(dccnpath('/home/common/matlab/fieldtrip/data/test/latest/raw/meg/preproc_ctf151.mat'));

cfg = [];
cfg.method = 'summary';
cfg.layout = 'CTF151.lay';
data = ft_rejectvisual(cfg, data);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
2*((x - x2)*((x - x1)^2 + (y - y1)^2)*(-r12^2 + (r11 -
    ((x - x1)^2 + (y - y1)^2)^(1/2))^2 + (z - z1)^2)*(-r22^2 + (y -
    y2)^2 + (r21 - ((x - x2)^2 + (z - z2)^2)^(1/2))^2) + (x -
    x2)*((x - x1)^2 + (y - y1)^2)*(-r12^2 + (r11 -
    ((x - x1)^2 + (y - y1)^2)^(1/2))^2 + (z - z1)^2)*(-r21 +
    ((x - x2)^2 + (z - z2)^2)^(1/2))*((x - x2)^2 + (z -
    z2)^2)^(1/2) + (x - x1)*(-r11 +
    ((x - x1)^2 + (y - y1)^2)^(1/2))*((x - x1)^2 + (y -
    y1)^2)^(1/2)*(-r22^2 + (y - y2)^2 + (r21 -
    ((x - x2)^2 + (z - z2)^2)^(1/2))^2)*((x - x2)^2 + (z -
    z2)^2) + (x - x1)*(-r12^2 + (r11 - ((x - x1)^2 + (y - y1)^2)^(1/2))^2 +
    (z - z1)^2)*(-r22^2 + (y - y2)^2 + (r21 -
    ((x - x2)^2 + (z - z2)^2)^(1/2))^(2))*((x - x2)^2 + (z - z2)^2));

2*(((x - x1)^2 + (y - y1)^2)*(y -
    y2)*(-r12^2 + (r11 - ((x - x1)^2 + (y - y1)^2)^(1/2))^2 +
    (z - z1)^2) + (-r11 + ((x - x1)^2 + (y - y1)^2)^(1/2))*
    ((x - x1)^2 + (y - y1)^2)^(1/2)*(y - y1)*(-r22^2 +
    (y - y2)^2 + (r21 - ((x - x2)^2 + (z - z2)^2)^(1/2))^2) +
    (y - y1)*(-r12^2 + (r11 - ((x - x1)^2 + (y - y1)^2)^(1/2))^2 +
    (z - z1)^2)*(-r22^2 + (y - y2)^2 + (r21 -
    ((x - x2)^2 + (z - z2)^2)^(1/2))^2))*((x - x2)^2 + (z - z2)^2);

2*((x - x1)^2 + (y - y1)^2)*((z - z1)*(-r22^2 + (y - y2)^2 +
    (r21 - ((x - x2)^2 + (z - z2)^2)^(1/2))^2)*((x - x2)^2 +
    (z - z2)^2) + (-r12^2 + (r11 - ((x - x1)^2 + (y - y1)^2)^(1/2))^2 +
    (z - z1)^2)*(-r22^2 + (y - y2)^2 + (r21 - ((x - x2)^2 + (z - z2)^2)^(1/2))^2)*
    (z - z2) + (-r12^2 + (r11 - ((x - x1)^2 + (y - y1)^2)^(1/2))^2+
    (z - z1)^2)*(-r21 + ((x - x2)^2 + (z - z2)^2)^(1/2))*
    ((x - x2)^2 + (z - z2)^2)^(1/2)*(z - z2));

repo: ifem


text:
function [ t, dtdr, dtds ] = shape_t10 ( r, s )

%*****************************************************************************80
%
%% SHAPE_T10 evaluates shape functions for a 10 node reference triangle.
%
%  Reference Element T10:
%
%    |
%    1  10
%    |  |\
%    |  | \
%    |  8  9
%    |  |   \
%    S  |    \
%    |  5  6  7
%    |  |      \
%    |  |       \
%    0  1--2--3--4
%    |
%    +--0----R---1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, S, the reference coordinates of a point.
%
%    Output, real T(10), the basis functions at the point.
%
%    Output, real DTDR(10), the R basis derivatives at the point.
%
%    Output, real DTDS(10), the S basis derivatives at the point.
%
  a = 1.0 / 3.0;
  b = 2.0 / 3.0;
  c = 1.0;

  t(1)  = 4.5 * ( a - r - s ) * ( b - r - s ) * ( c - r - s );
  t(2)  = 13.5 * r * ( b - r - s ) * ( c - r - s );
  t(3)  = - 13.5 * r * ( a - r ) * ( c - r - s );
  t(4)  = 4.5 * r * ( a - r ) * ( b - r );
  t(5)  = 13.5 * s * ( b - r - s ) * ( c - r - s );
  t(6)  = 27.0 * r * s * ( c - r - s );
  t(7)  = 13.5 * r * s * ( r - a );
  t(8)  = 13.5 * s * ( s - a ) * ( c - r - s );
  t(9)  = 13.5 * r * s * ( s - a );
  t(10) = 4.5 * s * ( a - s ) * ( b - s );

  dtdr(1) = 4.5 * ( ( a - s ) * ( 2.0 * r - c - b + 2.0 * s ) ...
    - ( s - b ) * ( s - c ) - 2.0 * ( 2.0 * s - b - c ) * r - 3.0 * r * r );
  dtdr(2) = 13.5 * ( ( s - b ) * ( s - c ) + 2.0 * ( 2.0 * s - b - c ) * r ...
    + 3.0 * r * r );
  dtdr(3) = - 13.5 * ( a * ( c - s ) + 2.0 * ( s - a - c ) * r + 3.0 * r * r );
  dtdr(4) = 4.5 * ( a * b - 2.0 * ( a + b ) * r + 3.0 * r * r );
  dtdr(5) = 13.5 * s * ( 2.0 * s - b - c + 2.0 * r );
  dtdr(6) = 27.0 * s * ( c - s - 2.0 * r );
  dtdr(7) = 13.5 * s * ( 2.0 * r - a );
  dtdr(8) = - 13.5 * s * ( s - a );
  dtdr(9) = 13.5 * s * ( s - a);
  dtdr(10) = 0.0;

  dtds(1) = 4.5 * ( ( a - r ) * ( 2.0 * s - c - b + 2.0 * r ) ...
    - ( r - b ) * ( r - c ) - 2.0 * ( 2.0 * r - b - c ) * s - 3.0 * s * s );
  dtds(2) = 13.5 * r * ( 2.0 * s + 2.0 * r - b - c );
  dtds(3) = 13.5 * r * ( a - r );
  dtds(4) = 0.0;
  dtds(5) = 13.5 * ( ( r - b ) * ( r - c ) + 2.0 * ( 2.0 * r - b - c ) * s ...
    + 3.0 * s * s );
  dtds(6) = 27.0 * r * ( c - r - 2.0 * s );
  dtds(7) = 13.5 * r * ( r - a );
  dtds(8) = - 13.5 * ( a * ( c - r ) + 2.0 * ( r - c - a ) * s + 3.0 * s * s );
  dtds(9) = 13.5 * r * ( 2.0 * s - a );
  dtds(10) = 4.5 * ( a * b - 2.0 * ( a + b ) * s + 3.0 * s * s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function boundary = component_boundary(grid, ndx)

if length(ndx) == 2
  ndx = sub2ind(size(grid), ndx(1), ndx(2));
end

filled = imfill(~grid, ndx, 8);
component = ~(xor(filled,grid));
boundary = iris.terrain_grid.find_boundary(component);

% figure(21)
% subplot(411)
% imshow(grid)
% subplot(412)
% imshow(component)
% subplot(413)
% imshow(fast_boundary);
% subplot(414)
% imshow(boundary);
% assert(all(all(fast_boundary == boundary)));

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
function model = yalmip2dsdp(interfacedata);

[C,A,b,blk] = sedumi2dsdp(interfacedata.F_struc,interfacedata.c,interfacedata.K);
interfacedata.options.dsdp.dual_quadratic=spalloc(length(interfacedata.c),length(interfacedata.c),0);
interfacedata.options.dsdp.printyes = (interfacedata.options.verbose>0);
model.A = A;
model.C = C;
model.b = b;
model.options = interfacedata.options.dsdp

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function this = save(this)
%saves MrSeries in different file formats
%
%   MrSeries = save(MrSeries)
%
% This is a method of class MrSeries.
%
% IN
%   parameters.save.format
%
% OUT
%
% EXAMPLE
%   save
%
%   See also MrSeries

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-07-02
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

pathSave = this.parameters.save.path;

if ~exist(pathSave, 'dir')
    mkdir(pathSave)
end

% set save path of all objects to the one of the MrSeries
this.set_save_path();

% detect all image files from object and save as nifti-files in save-path
% of MrSeries
handleImageArray = this.get_all_image_objects();
for iImage = 1:numel(handleImageArray);
    handleImageArray{iImage}.save;
end

% strip data from object and save MrSeries itself
MrObject = this.copyobj('exclude', 'data'); % copies object without data
fileObject = fullfile(pathSave, 'MrObject.mat');
save(fileObject, 'MrObject')

author: translationalneuromodeling

repo: tapas


text:
function [x,mc,mn,mx]=v_melbankm(p,n,fs,fl,fh,w)
%V_MELBANKM determine matrix for a mel/erb/bark-spaced filterbank [X,MN,MX]=(P,N,FS,FL,FH,W)
%
% Inputs:
%       p   number of filters in v_filterbank or the filter spacing in k-mel/bark/erb [ceil(4.6*log10(fs))]
%		n   length of fft
%		fs  sample rate in Hz
%		fl  low end of the lowest filter as a fraction of fs [default = 0]
%		fh  high end of highest filter as a fraction of fs [default = 0.5]
%		w   any sensible combination of the following:
%             'b' = bark scale instead of mel
%             'e' = erb-rate scale
%             'l' = log10 Hz frequency scale
%             'f' = linear frequency scale
%
%             'c' = fl/fh specify centre of low and high filters
%             'h' = fl/fh are in Hz instead of fractions of fs
%             'H' = fl/fh are in mel/erb/bark/log10
%
%		      't' = triangular shaped filters in mel/erb/bark domain (default)
%		      'n' = hanning shaped filters in mel/erb/bark domain
%		      'm' = hamming shaped filters in mel/erb/bark domain
%
%		      'z' = highest and lowest filters taper down to zero [default]
%		      'y' = lowest filter remains at 1 down to 0 frequency and highest filter
%                   remains at 1 up to nyquist freqency. See note (1) below.
%
%             'u' = scale filters to sum to unity
%
%             's' = single-sided: do not double filters to account for negative frequencies
%
%             'g' = plot idealized filters [default if no output arguments present]
%
% Outputs:	x     a sparse matrix containing the v_filterbank amplitudes
%		          If the mn and mx outputs are given then size(x)=[p,mx-mn+1]
%                 otherwise size(x)=[p,1+floor(n/2)]
%                 Note that the peak filter values equal 2 to account for the power in the negative FFT frequencies.
%           mc    the v_filterbank centre frequencies in mel/erb/bark
%		    mn    the lowest fft bin with a non-zero coefficient
%		    mx    the highest fft bin with a non-zero coefficient
%                 NOTE: For legacy compatibility reasons, you must specify both or neither of mn and mx.
%
% Notes: (1) If 'ty' or 'ny' is specified, the total power in the fft is preserved.
%        (2) The filter shape (triangular, hamming etc) is defined in the mel (or erb etc) domain
%            rather than in the linear frequency domain which is more common (e.g. [2]).
%        (3) A mel-filterbank can also be created using v_filtbank() which uses triangular
%            filters in the linear frequency domain and copes better with the narrow filters
%            that arise when p is large on n is small.
%
% Examples of use:
%
% (a) Calcuate the Mel-frequency Cepstral Coefficients
%
%       f=v_rfft(s);			        % v_rfft() returns only 1+floor(n/2) coefficients
%		x=v_melbankm(p,n,fs);	        % n is the fft length, p is the number of filters wanted
%		z=log(x*abs(f).^2);         % multiply x by the power spectrum
%		c=dct(z);                   % take the DCT
%
% (b) Calcuate the Mel-frequency Cepstral Coefficients efficiently
%
%       f=fft(s);                        % n is the fft length, p is the number of filters wanted
%       [x,mc,na,nb]=v_melbankm(p,n,fs);   % na:nb gives the fft bins that are needed
%       z=log(x*(f(na:nb)).*conj(f(na:nb)));
%		c=dct(z);                   % take the DCT
%
% (c) Plot the calculated filterbanks
%
%      plot((0:floor(n/2))*fs/n,melbankm(p,n,fs)')   % fs=sample frequency
%
% (d) Plot the idealized filterbanks (without output sampling)
%
%      v_melbankm(p,n,fs);
%
% References:
%
% [1] S. S. Stevens, J. Volkman, and E. B. Newman. A scale for the measurement
%     of the psychological magnitude of pitch. J. Acoust Soc Amer, 8: 185-19, 1937.
% [2] S. Davis and P. Mermelstein. Comparison of parametric representations for
%     monosyllabic word recognition in continuously spoken sentences.
%     IEEE Trans Acoustics Speech and Signal Processing, 28 (4): 357-366, Aug. 1980.

%      Copyright (C) Mike Brookes 1997-2009
%      Version: $Id: v_melbankm.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Note, in the comments, "FFT bin_0" assumes DC = bin 0 whereas "FFT bin_1" assumes DC = bin 1

if nargin < 6
    w='tz'; % default options
end
if nargin < 5 || isempty(fh)
    fh=0.5; % max freq is the nyquist
end
if nargin < 4 || isempty(fl)
    fl=0; % min freq is DC
end

sfact=2-any(w=='s');   % 1 if single sided else 2
wr=' ';   % default warping is mel
for i=1:length(w)
    if any(w(i)=='lebf')
        wr=w(i);
    end
end
if any(w=='h') || any(w=='H')
    mflh=[fl fh];
else
    mflh=[fl fh]*fs;
end
if ~any(w=='H')
    switch wr
        case 'f'                        % no transformation
        case 'l'
            if fl<=0
                error('Low frequency limit must be >0 for ''l'' option');
            end
            mflh=log10(mflh);           % convert frequency limits into log10 Hz
        case 'e'
            mflh=v_frq2erb(mflh);       % convert frequency limits into erb-rate
        case 'b'
            mflh=v_frq2bark(mflh);      % convert frequency limits into bark
        otherwise
            mflh=v_frq2mel(mflh);       % convert frequency limits into mel
    end
end
melrng=mflh*(-1:2:1)';                  % mel range
fn2=floor(n/2);                         % bin index of highest positive frequency (Nyquist if n is even)
if isempty(p)
    p=ceil(4.6*log10(fs));              % default number of filters
end
if any(w=='c')                          % c option: specify filter centres not edges
    if p<1
        p=round(melrng/(p*1000))+1;
    end
    melinc=melrng/(p-1);
    mflh=mflh+(-1:2:1)*melinc;
else
    if p<1
        p=round(melrng/(p*1000))-1;
    end
    melinc=melrng/(p+1);
end
%
% Calculate the FFT bins corresponding to [filter#1-low filter#1-mid filter#p-mid filter#p-high]
%
switch wr
    case 'f'
        blim=(mflh(1)+[0 1 p p+1]*melinc)*n/fs;
    case 'l'
        blim=10.^(mflh(1)+[0 1 p p+1]*melinc)*n/fs;
    case 'e'
        blim=v_erb2frq(mflh(1)+[0 1 p p+1]*melinc)*n/fs;
    case 'b'
        blim=v_bark2frq(mflh(1)+[0 1 p p+1]*melinc)*n/fs;
    otherwise
        blim=v_mel2frq(mflh(1)+[0 1 p p+1]*melinc)*n/fs;
end
mc=mflh(1)+(1:p)*melinc;        % mel centre frequencies
b1=floor(blim(1))+1;            % lowest FFT bin_0 required (might be negative)
b4=min(fn2,ceil(blim(4))-1);    % highest FFT bin_0 required
%
% now map all the useful FFT bins_0 to filter_1 centres
%
switch wr
    case 'f'
        pf=((b1:b4)*fs/n-mflh(1))/melinc;
    case 'l'
        pf=(log10((b1:b4)*fs/n)-mflh(1))/melinc;
    case 'e'
        pf=(v_frq2erb((b1:b4)*fs/n)-mflh(1))/melinc;
    case 'b'
        pf=(v_frq2bark((b1:b4)*fs/n)-mflh(1))/melinc;
    otherwise
        pf=(v_frq2mel((b1:b4)*fs/n)-mflh(1))/melinc;
end
%
%  remove any incorrect entries in pf due to rounding errors
%
if pf(1)<0
    pf(1)=[];
    b1=b1+1;
end
if pf(end)>=p+1
    pf(end)=[];
    b4=b4-1;
end
fp=floor(pf);                   % FFT bin_0 i contributes to filters_1 fp(1+i-b1)+[0 1]
pm=pf-fp;                       % multiplier for upper filter
k2=find(fp>0,1);                % FFT bin_1 k2+b1 is the first to contribute to both upper and lower filters
k3=find(fp<p,1,'last');         % FFT bin_1 k3+b1 is the last to contribute to both upper and lower filters
k4=numel(fp);                   % FFT bin_1 k4+b1 is the last to contribute to any filters
if isempty(k2)
    k2=k4+1;
end
if isempty(k3)
    k3=0;
end
if any(w=='y')                  % preserve power in FFT
    mn=1;                       % lowest fft bin required (1 = DC)
    mx=fn2+1;                   % highest fft bin required (1 = DC)
    r=[ones(1,k2+b1-1) 1+fp(k2:k3) fp(k2:k3) repmat(p,1,fn2-k3-b1+1)];  % filter number_1
    c=[1:k2+b1-1 k2+b1:k3+b1 k2+b1:k3+b1 k3+b1+1:fn2+1];                % FFT bin1
    v=[ones(1,k2+b1-1) pm(k2:k3) 1-pm(k2:k3) ones(1,fn2-k3-b1+1)];
else
    r=[1+fp(1:k3) fp(k2:k4)];       % filter number_1
    c=[1:k3 k2:k4];                 % FFT bin_1 - b1
    v=[pm(1:k3) 1-pm(k2:k4)];
    mn=b1+1;                        % lowest fft bin_1
    mx=b4+1;                        % highest fft bin_1
end
if b1<0
    c=abs(c+b1-1)-b1+1;             % convert negative frequencies into positive
end
% end
if any(w=='n')
    v=0.5-0.5*cos(v*pi);            % convert triangles to Hanning
elseif any(w=='m')
    v=0.5-0.46/1.08*cos(v*pi);      % convert triangles to Hamming
end
if sfact==2                         % double all except the DC and Nyquist (if any) terms
    msk=(c+mn>2) & (c+mn<n-fn2+2);  % there is no Nyquist term if n is odd
    v(msk)=2*v(msk);
end
%
% sort out the output argument options
%
if nargout > 2
    x=sparse(r,c,v);
    if nargout == 3                 % if exactly three output arguments, then
        mc=mn;                      % delete mc output for legacy code compatibility
        mn=mx;
    end
else
    x=sparse(r,c+mn-1,v,p,1+fn2);
end
if any(w=='u')
    sx=sum(x,2);
    x=x./repmat(sx+(sx==0),1,size(x,2));
end
%
% plot results if no output arguments or g option given
%
if ~nargout || any(w=='g')          % plot idealized filters
    ng=201;                         % 201 points
    me=mflh(1)+(0:p+1)'*melinc;
    switch wr
        case 'f'
            fe=me;                  % defining frequencies
            xg=repmat(linspace(0,1,ng),p,1).*repmat(me(3:end)-me(1:end-2),1,ng)+repmat(me(1:end-2),1,ng);
        case 'l'
            fe=10.^me;              % defining frequencies
            xg=10.^(repmat(linspace(0,1,ng),p,1).*repmat(me(3:end)-me(1:end-2),1,ng)+repmat(me(1:end-2),1,ng));
        case 'e'
            fe=v_erb2frq(me);       % defining frequencies
            xg=v_erb2frq(repmat(linspace(0,1,ng),p,1).*repmat(me(3:end)-me(1:end-2),1,ng)+repmat(me(1:end-2),1,ng));
        case 'b'
            fe=v_bark2frq(me);      % defining frequencies
            xg=v_bark2frq(repmat(linspace(0,1,ng),p,1).*repmat(me(3:end)-me(1:end-2),1,ng)+repmat(me(1:end-2),1,ng));
        otherwise
            fe=v_mel2frq(me);       % defining frequencies
            xg=v_mel2frq(repmat(linspace(0,1,ng),p,1).*repmat(me(3:end)-me(1:end-2),1,ng)+repmat(me(1:end-2),1,ng));
    end

    v=1-abs(linspace(-1,1,ng));
    if any(w=='n')
        v=0.5-0.5*cos(v*pi);        % convert triangles to Hanning
    elseif any(w=='m')
        v=0.5-0.46/1.08*cos(v*pi);  % convert triangles to Hamming
    end
    v=v*sfact;                      % multiply by 2 if double sided
    v=repmat(v,p,1);
    if any(w=='y')                  % extend first and last filters
        v(1,xg(1,:)<fe(2))=sfact;
        v(end,xg(end,:)>fe(p+1))=sfact;
    end
    if any(w=='u')                  % scale to unity sum
        dx=(xg(:,3:end)-xg(:,1:end-2))/2;
        dx=dx(:,[1 1:ng-2 ng-2]);
        vs=sum(v.*dx,2);
        v=v./repmat(vs+(vs==0),1,ng)*fs/n;
    end
    plot(xg',v','b');
    set(gca,'xlim',[fe(1) fe(end)]);
    xlabel(['Frequency (' v_xticksi 'Hz)']);
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%  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
%TFDEMO4 Cohen's class time-frequency distributions.
%	Time-Frequency Toolbox demonstration.
%
%	See also TFDEMO.

%	O. Lemoine - May 1996.
%	Copyright (c) CNRS.

clc; zoom on; clf;
echo on;

% The Wigner-Ville distribution
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% A time-frequency energy distribution which is particularly interesting
% is the Wigner-Ville distribution (WVD),which satisfies a large number
% of desirable mathematical properties. Let us see what we obtain on two
% particular synthetic signals :
% - the first signal is the academic linear chirp signal :

sig=fmlin(128);

% If we choose a 3-dimension plot to represent it, we can see that the WVD
% can take negative values, and that the localization obtained in the
% time-frequency plane for this signal is almost perfect.

tfrwv(sig);
% Press any key to continue...

pause ; clc;

% - the second one illustrates the Doppler effect, which expresses the
% dependence of the frequency received by an observer from a transmitter
% on the relative speed between the observer and the transmitter :

[fm,am,iflaw]=doppler(128,50,13,10,200);
sig=am.*fm;
tfrwv(sig);
% Looking at this time-frequency distribution, we notice that the energy is
% not distributed as we could expect for this signal. Although the signal
% term is well localized in the time-frequency plane, numerous other terms
% (the interference terms, due to the bilinearity of the WVD) are present at
% positions in time and frequency where the energy should be null.
%
% Press any key to continue...

pause; clc; close

% Interference geometry of the WVD
%""""""""""""""""""""""""""""""""""
% The rule of interference construction of the WVD can be summarized as
% follows : two points of the time-frequency plane interfere to create a
% contribution on a third point which is located at their geometrical
% midpoint. Besides, these interference terms oscillate perpendicularly to
% the line joining the two points interfering, with a frequency
% proportional to the distance between these two points.
%  This can be seen on the following example : we consider two atoms in
% the time-frequency plane, analyzed by the WVD, whose relative distance
% is increasing from one realization to the other, and then decreasing.
% The WVDs were calculated and saved on the file movieat.mat. We load them
% and run the sequence using the function movie :

M=movwv2at(128,15); movie(M,10); clear M;

% We can notice, from this movie, the evolution of the interferences
% when the distance between the two interfering terms changes, and in
% particular the change in the period and the direction of the oscillations.
%
% Press any key to continue...

pause; clc;

% The pseudo-WVD
%""""""""""""""""
%   As the analyzed signal is not known from -infinity to +infinity in
% pratice, one often consider a windowed version of the WVD, called the
% pseudo-WVD. The time-windowing operated has the effect of smoothing the
% WVD in frequency. Thus, because of their oscillating nature, the
% interferences will be attenuated in the pseudo-WVD compared to the WVD.
% However, the consequence of this improved readability is that many
% properties of the WVD are lost.
% If we consider a signal composed of four gaussian atoms, each localized
% at a corner of a rectangle,

sig=atoms(128,[32,.15,20,1;96,.15,20,1;32,.35,20,1;96,.35,20,1]);

% and compute its WVD

tfrwv(sig);
% we can see the four signal terms, along with six interference terms (two of
% them are superimposed). If we now compute the pseudo-WVD,

figure
tfrpwv(sig);
% we can note the important attenuation of the interferences oscillating
% perpendicularly to the frequency axis, and in return the spreading in
% frequency of the signal terms.
%
% Press any key to continue...

pause; clc; close;

% Importance of the analytic signal
%"""""""""""""""""""""""""""""""""""
%   Due to the quadratic nature of the WVD, its discrete version may be
% affected by a spectral aliasing, in particular if the signal x is
% real-valued and sampled at the Nyquist rate. A solution to this problem
% consists in using the analytic signal. Indeed, as its bandwidth is half the
% one of the real signal, the aliasing will not take place in the useful
% spectral domain [0,1/2] of this signal. This solution presents a second
% advantage : since the spectral domain is divided by two, the number of
% components in the time-frequency plane is also divided by two.
% Consequently, the number of interference terms decreases significantly.
% Here is an illustration : we first consider the WVD of the real part of
% a signal composed of two atoms :

sig=atoms(128,[32,0.15,20,1;96,0.32,20,1]);
tfrwv(real(sig));
% We can see that four signal terms are present instead of two, due to the
% spectral aliasing. Besides, because of the components located at negative
% frequencies (between -1/2 and 0), additional interference terms are
% present.
%
% Press any key to continue...

pause;

% If we now consider the WVD of the same signal, but in its complex
% analytic form,

tfrwv(sig);
% the aliasing effect has disappeared, as well as the terms corresponding to
% interferences between negative- and positive- frequency components.
%
% Press any key to continue...

pause; clc;

% The Cohen's class
%~~~~~~~~~~~~~~~~~~~
%  The Cohen's class gather all the time-frequency energy distributions which
% are covariant by translations in time and in frequency. It can be expressed
% as a 2-D correlation between a function PI(t,nu) and the WVD of the
% analyzed signal. The WVD is the element of the Cohen's class for which PI
% is a double Dirac, and the spectrogram is th element for which PI is the
% WVD of the short time window h. We consider in the following other elements
% of this class
%
% The smoothed pseudo-WVD
%"""""""""""""""""""""""""
% If we consider a separable smoothing function PI(t,nu)=g(t)H(-nu) (where
% H(nu) is the Fourier transform of a smoothing window h(t)), we allow a
% progressive and independent control, in both time and frequency, of the
% smoothing applied to the WVD. The obtained distribution is known as the
% smoothed-pseudo WVD.
% For example,let's consider a signal composed of two components : the first
% one is a complex sinusoid (normalized frequency 0.15) and the second one
% is a Gaussian signal shifted in time and frequency :

sig=sigmerge(fmconst(128,.15),amgauss(128).*fmconst(128,0.4),5);

% If we display the WVD, the pseudo-WV and the smoothed-pseudo-WVD of it,

tfrwv(sig);
figure; tfrpwv(sig);
figure; tfrspwv(sig);
% we can make the following remarks : from the WVD, we can see the two
% signal terms located at the right positions in the time-frequency plane,
% as well as the interference terms between them. As these interference
% terms oscillate globally perpendicularly to the time-axis, the frequency
% smoothing done by the pseudo-WVD degrades the frequency resolution without
% really attenuating the interferences. On the other hand, the time-smoothing
% carried out by the smoothed-pseudo-WVD considerably reduces these
% interferences.
%
% Press any key to continue...

pause; clc; close; close; close

%   An interresting property of the smoothed-pseudo WVD is that it allows a
% continuous passage from the spectrogram to the WVD, under the condition
% that the smoothing functions g and h are gaussian. This is clearly
% illustrated by the function movsp2wv.m, which considers different
% transitions, on a signal composed of four atoms :

M=movsp2wv(128,15); movie(M,10); clear M;

% This movie shows the effect of a (time/frequency) smoothing on the
% interferences and on the resolutions : the WVD gives the best resolutions
% (in time and in frequency), but presents the most important interferences,
% whereas the spectrogram gives the worst resolutions, but with nearly no
% interferences ; and the smoothed-pseudo WVD allows to choose the best
% compromise between these two extremes.
%
% Press any key to continue...

pause; clc;

% The narrow-band ambiguity function
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  The narrow-band ambiguity function, often used in radar signal processing,
% is the two-dimensional Fourier transform of the WVD. This property can be
% used to attenuate some of the interference terms. Indeed, in the case of
% a multi-component signal, the elements of the AF corresponding to the
% signal components (denoted as the AF-signal terms) are mainly located
% around the origin, whereas the elements corresponding to interferences
% between the signal components (AF-interference terms) appear at a distance
% from the origin which is proportional to the time-frequency distance
% between the involved components. This can be noticed on a simple example :
% We apply consider a signal composed of two linear FM signals with
% gaussian amplitudes :

N=64; sig1=fmlin(N,0.2,0.5).*amgauss(N);
sig2=fmlin(N,0.3,0).*amgauss(N);
sig=[sig1;sig2];

% Let us first have a look at the WVD :

tfrwv(sig);
% We have two distinct signal terms, and some interferences oscillating in
% the middle.
%
% Press any key to continue...

pause;

% If we look at the ambiguity function of this signal,

clf; ambifunb(sig);

% we have around the origin (in the middle of the image) the AF-signal
% terms, whereas the AF-interference terms are located away from the origin.
%
% Press any key to continue...

pause; clc;

% Other important energy distributions
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% The Rihaczek and Margenau-Hill distributions
% """""""""""""""""""""""""""""""""""""""""""""
%  The Rihaczek distribution, defined as
%	Rx(t,nu)=x(t)X*(nu)e^{-j2pi nu t},
% is a complex energy density at point (t,nu). This distribution, which
% corresponds to the element of the Cohen's class for which
% f(xi,tau)=e^{jpi xi tau}, verifies many good properties. However, it is
% complex valued, which can be awkward in practice. The real part of the
% Rihaczek distribution is also a time-frequency distribution of the
% Cohen's class (f(xi,tau)=cos(pi xi tau)), known as the Margenau-Hill
% distribution
%   The interference structure of the Rihaczek and Margenau-Hill
% distributions is different from the Wigner-Ville one : the interference
% terms corresponding to two points located on (t1,nu1) and (t2,nu2) are
% positioned at the coordinates (t1,nu2) and (t2,nu1). This can be seen on
% the following example :

sig=atoms(128,[32,0.15,20,1;96,0.32,20,1]);
tfrmh(sig);
% Thus, the use of the Rihaczek (or Margenau-Hill) distribution for signals
% composed of multi-components located at the same position in time or in
% frequency is no advised, since the interference terms will then be
% superposed to the signal terms.
%
% Press any key to continue...

pause; clc; close

% The Choi-Williams distribution
% """""""""""""""""""""""""""""""
%  An example of reduced interference distribution is given by the
% Choi-Williams distribution, defined as
% CWx(t,nu)=sqrt(2/pi)\int\int sigma/|tau| e^{-2sigma^2(s-t)^2/tau^2}
%	    x(s+tau/2)x*(s-tau/2) e^{-j2pi nu tau} ds dtau
% Note that when sigma->+infty, we obtain the WVD. Inversely, the smaller
% sigma, the better the reduction of the interferences.
%   The "cross"-shape of the parametrization function of the Choi-Williams
% distribution implies that the efficiency of this distribution strongly
% depends on the nature of the analyzed signal. For instance, if the signal
% is composed of synchronized components in time or in frequency, the
% Choi-Williams distribution will present strong interferences. This can be
% observed on the following example : we analyze four gaussian atoms
% positionned at the corners of a rectangle rotating around the center of
% the time-frequency plane :

M=movcw4at(128); movie(M,5); clear M;

% When the time/frequency supports of the atoms overlap, some
% AF-interference terms are not be completly attenuated (those present
% around the axes of the ambiguity plane), and the efficiency of the
% distribution is quite poor.
%
% Press any key to continue...

pause; clc;

% Comparison of the parametrization functions
% """"""""""""""""""""""""""""""""""""""""""""
%  To illustrate the differences between some of the presented
% distributions, we represent their weighting (parametrization) function in
% the ambiguity plane, along with the result obtained by applying them on a
% two-component signal embedded in white gaussian noise : the signal is the
% sum of two linear FM signals, the first one with a frequency going from
% 0.05 to 0.15, and the second one from 0.2 to 0.5. The signal to noise
% ratio is 10 dB.

% On the left-hand side of the figures, the parametrization functions are
% represented in a schematic way by the bold contour lines (the weighting
% functions are mainly non-zeros inside these lines), superimposed to the
% ambiguity function of the signal. The AF-signal terms are in the middle of
% the ambiguity plane, whereas the AF-interference terms are distant from the
% center. On the right-hand side, the corresponding time-frequency
% distributions are represented.

paramfun

% From these plots, we can conclude that the ambiguity plane is very
% enlightening with regard to interference reduction in the case of
% multicomponent signals. On this example, we notice that the
% smoothed-pseudo-WVD is a particularly convenient and versatile
% candidate. This is due to the fact that we can adapt independently the
% time-width and frequency-width of its weighting function. But in the
% general case, it is interesting to have several distributions at our
% disposal since each one is well adapted to a certain type of
% signal. Besides, for a given signal, as a result of the different
% interference geometries, these distributions offer complementary
% descriptions of this signal.
%
% Press any key to end this demonstration

pause;
echo off

close;

author: HeLiangHIT


text:
function       [nDCnlX,blk_arr,DC,par] = Image2PGs_II( im, par)
% record the non-local patch set and the index of each patch in
% of seed patches in image
% Pad noisy image to avoid Borader Issues
Paddedim = padarray(im,[par.win,par.win],'symmetric','both');
X = zeros(par.ps2ch, par.maxrcp, 'double');
k = 0;
for l = 1:par.ch
    for i = 1:par.ps
        for j = 1:par.ps
            k = k+1;
            blk = Paddedim(i:end-par.ps+i,j:end-par.ps+j,l);
            X(k,:) = blk(:)';
        end
    end
end
% index of each patch in Pad image
Index = (1:par.maxrcp);
Index = reshape(Index,par.maxrp,par.maxcp);
% record the indexs of patches similar to the seed patch
blk_arr = zeros(1, par.lenrc*par.nlsp,'double');
% Patch Group Means
DC = zeros(par.ps2ch,par.lenrc*par.nlsp,'double');
% non-local patch groups
nDCnlX = zeros(par.ps2ch,par.lenrc *par.nlsp,'double');
% record the distance of compared patches to the reference patches
Vdis_rp = zeros((2*par.win+1)^2, par.lenrc);
% record the index of compared patches to the reference patches
Vidx_rp = Vdis_rp;
% Compute the Integral Image
v = Paddedim(1:par.h+par.win,1:par.w+par.win);
k = 0;
for dy = -par.win:par.win
    for dx = -par.win:par.win
        k = k + 1;
        % initial distance matrix in each iteration
        Mdis_rp = Inf*ones(par.lenr, par.lenc, 'single' );
        % Decide shift type, tx = vx+dx; ty = vy+dy
        t = zeros(size(v));
        if dx == 0 && dy == 0
            Midx_rp = Index(par.r+par.win,par.c+par.win);
            Vidx_rp(k,:) = Midx_rp(:);
            continue;
        elseif dx <= 0 && dy <= 0
            t(-dx+1:end,-dy+1:end) = v(1:end+dx,1:end+dy);
            a = 1-floor(dx/par.step);
            b = par.lenr;
            c = 1-floor(dy/par.step);
            d = par.lenc;
        elseif dx <= 0 && dy > 0
            t(-dx+1:end,1:end-dy) = v(1:end+dx,dy+1:end);
            ddy = (dy>1)*(2+ceil((dy-2)/par.step)) + (0<dy && dy<=1)*dy;
            a = 1 - floor(dx/par.step);
            b = par.lenr;
            c = 1;
            d = par.lenc - ddy;
        elseif dx > 0 && dy <= 0
            t(1:end-dx,-dy+1:end) = v(dx+1:end,1:end+dy);
            ddx = (dx>1)*(2+ceil((dx-2)/par.step)) + (0<dx && dx<=1)*dx;
            a = 1;
            b = par.lenr - ddx;
            c = 1-floor(dy/par.step);
            d = par.lenc;
        elseif dx > 0 && dy > 0
            t(1:end-dx,1:end-dy) = v(dx+1:end,dy+1:end);
            ddx = (dx>1)*(2+ceil((dx-2)/par.step)) + (0<dx && dx<=1)*dx;
            ddy = (dy>1)*(2+ceil((dy-2)/par.step)) + (0<dy && dy<=1)*dy;
            a = 1;
            b = par.lenr - ddx;
            c = 1;
            d = par.lenc - ddy;
        end
        % Create sqaured difference image
        diff = (v-t).^2;
        % Construct integral image along rows
        Sd = cumsum(diff,1);
        % Construct integral image along columns
        Sd = cumsum(Sd,2);
        % Obtaine the Square difference for each reference patches
        SqDist = Sd(par.x,par.y) + Sd(par.x0,par.y0) - Sd(par.x,par.y0) - Sd(par.x0,par.y);
        % Obtain the real value Square difference for suitable reference patches
        Mdis_rp( a:b,c:d ) = SqDist(a:b,c:d);
        % corresponding distance to reference patches
        Vdis_rp(k,:) = Mdis_rp(:);
        % corresponding index of compared patch to reference patches
        Midx_rp = Index(par.r+par.win+dx,par.c+par.win+dy);
        Vidx_rp(k,:) = Midx_rp(:);
    end
end
[~,Ind] = sort(Vdis_rp);
for i = 1:size(Vdis_rp,2)
    indc = Vidx_rp( Ind( 1:par.nlsp,i ),i );
    blk_arr(1,(i-1)*par.nlsp+1:i*par.nlsp) = indc;
    temp = X( : , indc );
    DC(:,(i-1)*par.nlsp+1:i*par.nlsp) = repmat(mean(temp,2),[1 par.nlsp]);
    nDCnlX(:,(i-1)*par.nlsp+1:i*par.nlsp) = temp - DC(:,(i-1)*par.nlsp+1:i*par.nlsp);
end
blk_arr = par.maxr*(floor(blk_arr/par.maxrp)-par.win) + mod(blk_arr,par.maxrp) - par.win;

author: csjunxu


text:
function J = cell2Image(I)

    numChannels = numel(I);
    for channelIdx = 1:numChannels
        J(:,:,channelIdx) = I{channelIdx};
    end

author: thomas-koehler

repo: SupER


text:
function [lf, vol] = eeg_leadfield4(R, elc, vol)

% EEG_LEADFIELD4 electric leadfield for a dipole in 4 concentric spheres
%
% [lf] = eeg_leadfield4(R, elc, vol)
%
% with input arguments
%   R          position of the dipole
%   elc        position of the electrodes
% and vol being a structure with the elements
%   vol.r      radius of the 4 spheres
%   vol.cond   conductivity of the 4 spheres
%   vol.t      constant factors for series expansion (optional)
%
% The center of the spheres should be at the origin.
%
% This implementation is adapted from
%   Lutkenhoner, Habilschrift 1992.
% The original reference is
%  Cuffin BN, Cohen D. Comparison of the magnetoencephalogram and electroencephalogram. Electroencephalogr Clin Neurophysiol. 1979 Aug;47(2):132-46.
%
% See also EEG_LEADFIELD4_PREPARE for precomputing the constant factors,
% which can save time when multiple leadfield computations are done.

% Copyright (C) 2002, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% sort the spheres from the smallest to the largest
[vol.r, indx] = sort(vol.r);
vol.cond      = vol.cond(indx);

% use more convenient names for the radii and conductivities
r1 = vol.r(1); c1 = vol.cond(1);
r2 = vol.r(2); c2 = vol.cond(2);
r3 = vol.r(3); c3 = vol.cond(3);
r4 = vol.r(4); c4 = vol.cond(4);

% check whether the electrode ly on the sphere, allowing 0.5% tolerance
dist = sqrt(sum(elc.^2,2));
if any(abs(dist-r4)>r4*0.005)
  ft_warning('electrodes do not ly on sphere surface -> using projection')
end
elc = r4 * elc ./ [dist dist dist];

% check whether the dipole is inside the brain [disabled for EEGLAB]
% if sqrt(sum(R.^2))>=r1
%   ft_error('dipole is outside the brain compartment');
% end

% rotate everything so that the dipole is along the pos. z-axis
% only if the dipole is not in the origin or along the positive z-axis
if R(1)~=0 || R(2)~=0
  % compute the rotation matrix
  % the inverse rotation matrix is the transposed of this one
  val1 = norm(R);
  val2 = norm(R(1:2));
  rot(1,1) = R(1) * R(3) / (val1 * val2);
  rot(1,2) = R(2) * R(3) / (val1 * val2);
  rot(1,3) = -1.0 * val2 / val1;
  rot(2,1) = -1.0 * R(2) / val2;
  rot(2,2) =        R(1) / val2;
  rot(2,3) =                  0;
  rot(3,:) = R ./ val1;
  % rotate the electrodes
  elc = elc*rot';
elseif R(1)==0 && R(2)==0 && R(3)<0
  % dipole on negative z-axis, this case is very simple: reflect on xy-plane
  elc(:,3) = -elc(:,3);
else
  % dipole is on positive z-axis, nothing has to be done
end

% compute the constant factors for the sphere configuration if needed
if ~isfield(vol, 't')
  vol.t = eeg_leadfield4_prepare(vol);
end

Nchans = size(elc,1);
lf     = zeros(Nchans,3);
Nmax   = length(vol.t);
n      = 1:Nmax;
f      = norm(R)/r4;        % following cuffin1979
% c      = r2/r4;       % following cuffin1979
% d      = r3/r4;       % following cuffin1979

% this code is to cross-validate the lutkenhoner and cuffin implementations
% [lut_t, cuf_t] = eeg_leadfield4_prepare(vol);
% lut_c = (2*n+1).^4.*f.^(n-1) ./ (lut_t.*4*pi*c4*r4^2);
% cuf_c = (2*n+1).^4.*f.^(n-1) .*(c*d).^(2.*n+1) ./ (cuf_t.*4*pi*c4*r4^2);

% given a fixed volume conductor, these only need to be computed once for all electrodes
const = (2*n+1).^4.*f.^(n-1) ./ (vol.t.*4*pi*c4*r4^2);

for i=1:Nchans
  % convert the position of the electrodes to spherical coordinates
  [phi, el] = cart2sph(elc(i,1), elc(i,2), elc(i,3));

  % change from colatitude to latitude and compute the cosine
  cos_theta = cos(pi/2-el);

  % the series summation starts at zero
  s_x = 0;
  s_z = 0;

  for n=1:Nmax
    P0  = plgndr(n,0,cos_theta);        % zero'th order Legendre
    P1  = plgndr(n,1,cos_theta);        % first order Legendre
    s_x = s_x + const(n)*P1/n;          % s_y is identical
    s_z = s_z + const(n)*P0;
  end

  lf(i,1) = -cos(phi) * s_x;
  lf(i,2) = -sin(phi) * s_x;            % s_y is identical to s_x
  lf(i,3) = 1         * s_z;
end

% apply the inverse rotation to the leadfield matrix
if R(1)~=0 || R(2)~=0
  lf = lf*rot;
elseif R(1)==0 && R(2)==0 && R(3)<0
  lf(:,3) = -lf(:,3);
end


text:
function [G] = spm_morlet_conv(G,w,dt,wnum)
% temporal convolution of complex spectral responses with Morlet envelope
% FORMAT [G] = spm_morlet_conv(G,w,dt,wnum)
%
% G      - (t x w x n x n) cross spectral density
% w      - Frequencies (Hz)
% dt     - sampling interval (sec)
% wnum   - Wavelet number: default = 2  s.d. = wnum/(2*pi*w)
%
% G      - convolved cross spectral density
%__________________________________________________________________________
%
% This routine simply smooths a cross spectral response to emulate a
% wavelet transform.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_morlet_conv.m 6857 2016-08-19 15:17:06Z karl $

% setup and defaults
%--------------------------------------------------------------------------
if nargin < 4, wnum = 8; end
[nt,nw,ni,nj] = size(G);
iw            = 1:nw;

% get (non-stationary) convolution matrix for frequencies
%--------------------------------------------------------------------------
f     = (-nw:nw)';
H     = zeros(nw,nw);
for i = 1:nw
    s      = w(i)/wnum;
    h      = exp(-f.^2/(2*s^2));
    h      = h(iw + nw - i);
    H(:,i) = h/sum(h);
end

% convolution over frequencies
%--------------------------------------------------------------------------
for i = 1:ni
    for j = 1:nj
        G(:,:,i,j) = G(:,:,i,j)*H;
    end
end

% convolution over time
%--------------------------------------------------------------------------
for k = 1:nw
    s     = wnum/(2*pi*w(k));
    t     = -(s*4):dt:(s*4);
    h     = exp(-t.^2/(2*s^2));
    h     = spm_convmtx(h',nt,'square');
    h     = diag(1./sum(h,2))*h;
    for i = 1:ni
        for j = 1:nj
            G(:,k,i,j) = h*G(:,k,i,j);
        end
    end
end

author: spm


text:
function S = modexpon(S,p)
% Reduces the exponents of a sympoly modulo p
% Usage: S = modexpon(S,p)
%
% S is any sympoly
% p a scalar numeric value
%

if numel(S) > 1
  for i = 1:numel(S)
    S(i) = modexpon(S(i),p);
  end
else
  % a scalar sympoly
  S.Exponent = mod(S.Exponent,p);

  S = clean_sympoly(S);
end

author: Sable

repo: mcbench-benchmarks


text:
function x = blend_rs_1dn ( r, s, n, bound_rs )

%*****************************************************************************80
%
%% BLEND_RS_1DN extends vector data along sides into a square.
%
%  Diagram:
%
%    01-----r1-----11
%     |      .      |
%     |      .      |
%    0s.....rs.....1s
%     |      .      |
%     |      .      |
%    00-----r0-----10
%
%  Discussion:
%
%    BLEND_RS_1DN is NOT equivalent to a bilinear finite element method,
%    since the data is sampled everywhere along the boundary lines,
%    rather than at a finite number of nodes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Gordon,
%    Blending-Function Methods of Bivariate and Multivariate Interpolation
%    and Approximation,
%    SIAM Journal on Numerical Analysis,
%    Volume 8, Number 1, March 1971, pages 158-177.
%
%    William Gordon and Charles Hall,
%    Transfinite Element Methods: Blending-Function Interpolation over
%    Arbitrary Curved Element Domains,
%    Numerische Mathematik,
%    Volume 21, Number 1, 1973, pages 109-129.
%
%    William Gordon and Charles Hall,
%    Construction of Curvilinear Coordinate Systems and Application to
%    Mesh Generation,
%    International Journal of Numerical Methods in Engineering,
%    Volume 7, 1973, pages 461-477.
%
%    Joe Thompson, Bharat Soni, Nigel Weatherill,
%    Handbook of Grid Generation,
%    CRC Press, 1999.
%
%  Parameters:
%
%    Input, real R, S, the (R,S) coordinates of the point to be
%    evaluated.
%
%    Input, integer N, the dimension of the vector space.
%
%    External, BOUND_RS, is a function which is given (R,S)
%    coordinates and an component value I, and returns XI, the value
%    of the I-th component of the N-vector at that point.  BOUND_RS
%    will only be called for "sides", that is, for values (R,S) where
%    at least one of R and S is either 0.0 or 1.0.  BOUND_RS has the
%    form:
%      function xi = bound_rs ( r, s, i )
%
%    Output, real X(N), the interpolated value at the point (R,S).
%
  for i = 1 : n
%
%  Get the I-th coordinate component at the four corners.
%
    x00 = bound_rs ( 0.0, 0.0, i );
    x01 = bound_rs ( 0.0, 1.0, i );
    x10 = bound_rs ( 1.0, 0.0, i );
    x11 = bound_rs ( 1.0, 1.0, i );
%
%  Get the I-th coordinate component at the sides.
%
    xr0 = bound_rs ( r, 0.0, i );
    xr1 = bound_rs ( r, 1.0, i );
    x0s = bound_rs ( 0.0, s, i );
    x1s = bound_rs ( 1.0, s, i );
%
%  Interpolate the I-th coordinate component of the interior point.
%
    x(i) = blend_112 ( r, s, x00, x01, x10, x11, xr0, xr1, x0s, x1s );

  end

  return
end


text:
function [ n_data, n, m, x, fx ] = legendre_associated_normalized_values ...
  ( n_data )

%*****************************************************************************80
%
%% LEGENDRE_ASSOCIATED_NORMALIZED_VALUES: normalized associated Legendre.
%
%  Discussion:
%
%    The function considered is the associated Legendre polynomial P^M_N(X).
%
%    In Mathematica, the function can be evaluated by:
%
%      LegendreP [ n, m, x ]
%
%    The function is normalized by dividing by
%
%      sqrt ( 4 * pi * ( n + m )! / ( 2 * n + 1 ) / ( n - m )! )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Cambridge University Press, 1999,
%    ISBN: 0-521-64314-7,
%    LC: QA76.95.W65.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0
%    before the first call.  On each call, the routine increments N_DATA by 1,
%    and returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, integer M, real X,
%    the arguments of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 21;

  fx_vec = [ ...
     0.2820947917738781, ...
     0.2443012559514600, ...
    -0.2992067103010745, ...
    -0.07884789131313000, ...
    -0.3345232717786446, ...
     0.2897056515173922, ...
    -0.3265292910163510, ...
    -0.06997056236064664, ...
     0.3832445536624809, ...
    -0.2709948227475519, ...
    -0.2446290772414100, ...
     0.2560660384200185, ...
     0.1881693403754876, ...
    -0.4064922341213279, ...
     0.2489246395003027, ...
     0.08405804426339821, ...
     0.3293793022891428, ...
    -0.1588847984307093, ...
    -0.2808712959945307, ...
     0.4127948151484925, ...
    -0.2260970318780046 ]';
  m_vec = [ ...
    0, 0, 1, 0, ...
    1, 2, 0, 1, ...
    2, 3, 0, 1, ...
    2, 3, 4, 0, ...
    1, 2, 3, 4, ...
    5 ]';
  n_vec = [ ...
    0,  1,  1,  2, ...
    2,  2,  3,  3, ...
    3,  3,  4,  4, ...
    4,  4,  4,  5, ...
    5,  5,  5,  5, ...
    5 ]';
  x_vec = [ ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50, ...
    0.50 ]';

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    m = 0;
    x = 0.0;
    fx = 0.0;
  else
    n = n_vec(n_data);
    m = m_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function [Moda , ncont] = iModa(x)

% A matlab script to find the mode of vector x.
% There are two options at the end so that it is possible to perform a checksum
% on data and plot the vector of frequencies, just uncomment the right lines!
%
% Isaac M. M. INOGS, Trieste(Italy). January 10, 2005 @01h46:44
% Developed under Matlab(TM) version 7.0.0.19901 (R14)
% Modified on January 13, 2005 @14h49:50
% in order to include suggestion from Jos jnvdg@arcor.de (thanks, Jos)
% and add some robustness regarding NANs (Not-A-Number)

x = x(:);
Xvalue = unique(x(~isnan(x)));             % find elements in x with no repetitions (sorted output)

Xvalue = Xvalue(:);        % elements in vector
Xcont = histc(x,Xvalue);
Xcont = Xcont(:);          % frequency of each element

[maxCont imaxCont] = max(Xcont);
ncont = maxCont;                     % frequency of mode
Moda = Xvalue(imaxCont);             % the mode

%%%%%% optional 1 - plot of frequencies for each element

 plot(Xvalue,Xcont);
 xlabel('values')
 ylabel('frequencies')

%%%%%% optional 2 - checksum

% disp(strcat('sum of frequencies :',num2str(sum(Xcont))));
% disp(strcat('lenght of vector :',num2str(length(x))));
% disp(strcat('Nans, if any :',num2str(sum(isnan(x)))));
% chcksm = (sum(Xcont)+sum(isnan(x))) == length(x);
% disp(chcksm)


text:
function [valid, conditional] = valid_identifier(identifier)
% valid_identifier Verifies tracker identifier
%
% Tests if a given input string is a valid tracker identifier. A valid identifier contains only
% upper and lower case letters, digits and underscore.
%
% A conditionally valid identifier also contains dash. This is used for legacy reasons but is
% discouraged for new trackers.
%
% Input:
% - identifier: A string to test.
%
% Ouput:
% - valid: True if the string is a valid identifier.
% - conditional: True if the string is a conditionally valid identifier.

if ~ischar(identifier)
    valid = false;
	conditional = false;
    return;
end;

alpha = (identifier >= 'a' & identifier <= 'z') | (identifier >= 'A' & identifier <= 'Z');
numeric = (identifier >= '0' & identifier <= '9');
symbol = identifier == '_';

valid = all(alpha | numeric | symbol);
conditional = all(alpha | numeric | symbol | identifier == '-');

author: votchallenge

repo: toolkit-legacy


text:
function [X S Y dist] = OptSpace(M_E,r,niter,tol)
% An algorithm for Matrix Reconstruction from a partially revealed set.
% See "Matrix Completion from a Few Entries"(http://arxiv.org/pdf/0901.3150) for details
% Usage :
% [X S Y dist] = OptSpace(A,r,niter,tol);
% [X S Y dist] = OptSpace(A);
%
% INPUT :
% A     :  The partially revealed matrix.
%          Sparse matrix with zeroes at the unrevealed indices.
%
% r     :  The rank to be used for reconstruction. Use [] to guess the rank.
% niter :  The max. no. of iterations. Use [] to use default (50).
% tol   :  Stop iterations if norm( (XSY' - M_E).*E , 'fro' )/sqrt(|E|) < tol, where
%        - E_{ij} = 1 if M_{ij} is revealed and zero otherwise,
%        - |E| is the size of the revealed set.
%        - Use [] to use the default (1e-6)
%
%
% OUTPUT :
% X      : A size(A,1)xr matrix
% S      : An rxr matrix
% Y      : A size(A,2)xr matrix
% such that M_hat = X*S*Y'
% dist   : A vector containing norm( (XSY' - M_E).*E , 'fro' )/sqrt(|E|) at each
%          successive iteration
%
% Date : 21st April, 2009
% COPYRIGHT 2009 Raghunandan H. Keshavan, Andrea Montanari, Sewoong Oh

if(nargin==1)

	M_E = sparse(M_E);
	[n m] = size(M_E);
	E = spones(M_E);
	eps = nnz(E)/sqrt(m*n) ;

	tol = 1e-6;

	fprintf(1,'Rank not specified. Trying to guess ...\n');
	r = guessRank(M_E) ;
	fprintf(1,'Using Rank : %d\n',r);

	m0 = 10000 ;
	rho = 0;

	niter = 50;
elseif(nargin==4)

	M_E = sparse(M_E);
	[n m] = size(M_E);
	E = spones(M_E);
	eps = nnz(E)/sqrt(m*n) ;

	if( length(tol) == 0 )
		tol = 1e-6;
	end

	if( length(r) == 0 )

		fprintf(1,'Rank not specified. Trying to guess ...\n');
		r = guessRank(M_E) ;
		fprintf(1,'Using Rank : %d\n',r);
	end

	m0 = 10000 ;
	rho = 0;

	if( length(niter) == 0 )
		niter = 50 ;
	end
else
	fprintf(1,'Improper arguments (See "help OptSpace")\n');
	fprintf(1,'Usage :\n[X S Y dist] = OptSpace(A,r,niter,tol) \n') ;
	fprintf(1,'[X S Y dist] = OptSpace(A)\n');
	return;
end

rescal_param = sqrt( nnz(E) * r / norm(M_E,'fro')^2 ) ;
M_E = M_E * rescal_param ;

fprintf(1,'Trimming ...\n');
% Trimming

M_Et = M_E ;
d = sum(E);
d_=mean(full(d));
for col=1:m
    if ( sum(E(:,col))>2*d_ )
        list = find( E(:,col) > 0 );
        p = randperm(length(list));
        M_Et( list( p(ceil(2*d_):end) ) , col ) = 0;
    end
end

d = sum(E');
d_= mean(full(d));
for row=1:n
    if ( sum(E(row,:))>2*d_ )
        list = find( E(row,:) > 0 );
        p = randperm(length(list));
        M_Et(row,list( p(ceil(2*d_):end) ) ) = 0;
    end
end

fprintf(1,'Sparse SVD ...\n');
% Sparse SVD
[X0 S0 Y0] = svds(M_Et,r) ;

clear M_Et;

% Initial Guess
X0 = X0*sqrt(n) ; Y0 = Y0*sqrt(m) ;
S0 = S0 / eps ;

fprintf(1,'Iteration\tFit Error\n');

% Gradient Descent
X = X0;Y=Y0;
S = getoptS(X,Y,M_E,E);

dist(1) = norm( (M_E - X*S*Y').*E ,'fro')/sqrt(nnz(E) )  ;
fprintf(1,'0\t\t%e\n',dist(1) ) ;

for i = 1:niter

% Compute the Gradient
	[W Z] = gradF_t(X,Y,S,M_E,E,m0,rho);

% Line search for the optimum jump length
	t = getoptT(X,W,Y,Z,S,M_E,E,m0,rho) ;
	X = X + t*W;Y = Y + t*Z;S = getoptS(X,Y,M_E,E) ;

% Compute the distortion
	dist(i+1) = norm( (M_E - X*S*Y').*E,'fro' )/sqrt(nnz(E));
	fprintf(1,'%d\t\t%e\n',i,dist(i+1) ) ;
	if( dist(i+1) < tol )
		break ;
	end
end

S = S /rescal_param ;

% Function to Guess the Rank of the input Matrix
function r = guessRank(M_E);
	[n m] = size(M_E);
	epsilon = nnz(M_E)/sqrt(m*n);
    S0 = svds(M_E,100) ;

    S1=S0(1:end-1)-S0(2:end);
    S1_ = S1./mean(S1(end-10:end));
    r1=0;
    lam=0.05;
    while(r1<=0)
        for idx=1:length(S1_)
            cost(idx) = lam*max(S1_(idx:end)) + idx;
        end
        [v2 i2] = min(cost);
        r1 = max(i2-1);
        lam=lam+0.05;
    end

	clear cost;
    for idx=1:length(S0)-1
        cost(idx) = (S0(idx+1)+sqrt(idx*epsilon)*S0(1)/epsilon  )/S0(idx);
    end
    [v2 i2] = min(cost);
    r2 = max(i2);

	r = max([r1 r2]);

% * * * * * * * * * * * * * * * * * * * *
% Function to compute the distortion
function out = F_t(X,Y,S,M_E,E,m0,rho)
[n r] = size(X) ;

out1 = sum( sum( ( (X*S*Y' - M_E).*E ).^2 ) )/2 ;

out2 =  rho*G(Y,m0,r) ;
out3 =  rho*G(X,m0,r) ;
out = out1+out2+out3 ;

function out = G(X,m0,r)

z = sum(X.^2,2)/(2*m0*r) ;
y = exp( (z-1).^2 ) - 1 ;
y( find(z < 1) ) = 0 ;
out = sum(y) ;
% * * * * * * * * * * * * * * * * * * * *

% Function to compute the gradient
function [W Z] = gradF_t(X,Y,S,M_E,E,m0,rho)
[n r] = size(X);
[m r] = size(Y);

XS = X*S ;
YS = Y*S' ;
XSY = XS*Y' ;

Qx = X'* ( (M_E - XSY).*E )*YS /n;
Qy = Y'* ( (M_E - XSY).*E )'*XS /m;

W = ( (XSY - M_E).*E )*YS + X*Qx + rho*Gp(X,m0,r);
Z = ( (XSY - M_E).*E )'*XS + Y*Qy + rho*Gp(Y,m0,r);

function out = Gp(X,m0,r)
z = sum(X.^2,2) /(2*m0*r) ;
z = 2*exp( (z-1).^2 ).*(z-1) ;
z( find(z<0) ) = 0;

out = X.*repmat(z,1,r) / (m0*r) ;
% * * * * * * * * * * * * * * * * * * * *

% * * * * * * * * * * * * * * * * * * * *
% Function to find Sopt given X, Y
function out = getoptS(X,Y,M_E,E)

[n r] = size(X);
C = X' * ( M_E ) * Y ; C = C(:) ;

for i = 1:r
        for j = 1:r
                ind = (j-1)*r + i ;
                temp = X' * (  (X(:,i) * Y(:,j)').*E ) * Y ;
                A(:,ind) = temp(:) ;
        end
end

S = A\C ;
out = reshape(S,r,r) ;
% * * * * * * * * * * * * * * * * * * * *

% * * * * * * * * * * * * * * * * * * * *
% Function to perform line search
function out = getoptT(X,W,Y,Z,S,M_E,E,m0,rho)
norm2WZ = norm(W,'fro')^2 + norm(Z,'fro')^2;
f(1) = F_t(X, Y,S,M_E,E,m0,rho) ;

t = -1e-1 ;
for i = 1:20
        f(i+1) = F_t(X+t*W,Y+t*Z,S,M_E,E,m0,rho) ;

        if( f(i+1) - f(1) <= .5*(t)*norm2WZ )
            out = t ;
            return;
        end
        t = t/2 ;
end
out = t ;

author: andrewssobral

repo: lrslibrary


text:
function [DataOut, Metric] = ReceiverHD(HDchips, G, Scrambler);
%
% ReceiverHD    Hard Decision Receiver. This function performs transmitted packet data recovering,
%                          based on Viterbi Decoder.
%
% 						Block Diagram:
%						HDchips -> [DeScrambler] -> [De-Interleaver] -> [Viterbi Decoder]
%
%
% 						Inputs: HDchips - Hard Decision RAKE receiver symbols
%                         G - Viterbi Encoder generation polynom
%                         Scrambler - Scrambler sequence
%
% 						Outputs: DataOut - Received data (binary form)
%                                       Metric - The  best metric of Viterbi Decoder
%

%====================== COMMON DEFINITIONS =================
%------- Viterbi Polynom ---------
if (nargin == 1)
   G = [1 1 1 1 0 1 0 1 1; 1 0 1 1 1 0 0 0 1];
end

%============================== R E C E I V E R ==========================

%----------- DeScrambler---------
% Rate = 19.2 KBps
HDchips = xor(HDchips, Scrambler);

%-------- DeInterleaver ---------
INTERL = reshape(HDchips, 16, 24);		% IN-> rows, OUT-> columns
HDchips = reshape(INTERL', length(HDchips), 1);  % Rate = 19.2 KBps

%-------- HD Input Viterbi Decoder ----------
[DataOut Metric] = VitDec(G, HDchips, 1);


text:
function [chnsReg,chnsSim] = edgesChns( I, opts )
% Compute features for structured edge detection.
%
% For an introductory tutorial please see edgesDemo.m.
%
% USAGE
%  [chnsReg,chnsSim] = edgesChns( I, opts )
%
% INPUTS
%  I          - [h x w x 3] color input image
%  opts       - structured edge model options
%
% OUTPUTS
%  chnsReg    - [h x w x nChannel] regular output channels
%  chnsSim    - [h x w x nChannel] self-similarity output channels
%
% EXAMPLE
%
% See also edgesDemo, edgesTrain, edgesDetect, gradientMag
%
% Structured Edge Detection Toolbox      Version 3.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-microsoft.com]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the MSR-LA Full Rights License [see license.txt]

shrink=opts.shrink; nTypes=1; chns=cell(1,opts.nChns); k=0;
if(size(I,3)>3), nTypes=2; Is={I(:,:,1:3),I(:,:,4:end)}; end
for t=1:nTypes
  if(nTypes>1), I=Is{t}; end
  if(size(I,3)==1), cs='gray'; else cs='luv'; end; I=rgbConvert(I,cs);
  Ishrink=imResample(I,1/shrink); k=k+1; chns{k}=Ishrink;
  for i = 1:2, s=2^(i-1);
    if(s==shrink), I1=Ishrink; else I1=imResample(I,1/s); end
    I1 = convTri( I1, opts.grdSmooth );
    [M,O] = gradientMag( I1, 0, opts.normRad, .01 );
    H = gradientHist( M, O, max(1,shrink/s), opts.nOrients, 0 );
    k=k+1; chns{k}=imResample(M,s/shrink);
    k=k+1; chns{k}=imResample(H,max(1,s/shrink));
  end
end
chns=cat(3,chns{1:k}); assert(size(chns,3)==opts.nChns);
chnSm=opts.chnSmooth/shrink; if(chnSm>1), chnSm=round(chnSm); end
simSm=opts.simSmooth/shrink; if(simSm>1), simSm=round(simSm); end
chnsReg=convTri(chns,chnSm); chnsSim=convTri(chns,simSm);

end

author: Cloud-CV

repo: object-proposals


text:
function [xs,optimState] = searchCMA(sumrule,x,gpstruct,LB,UB,optimState,options)
%SEARCHCMA Covariance matrix adaptation search step (inspired by CMA-ES).

if nargin < 2 || isempty(x) || isempty(gpstruct)
    xs = 'cma';
    return;
end

% SUMRULE not specified, all arguments shifted by one
if nargin < 7
    options = optimState;
    optimState = UB;
    UB = LB;
    LB = gpstruct;
    gpstruct = x;
    x = sumrule;
    sumrule = [];
end

% By default renormalize by expected vector magnitude
if isempty(sumrule); sumrule = 1; end

MeshSize = optimState.meshsize;
SearchFactor = optimState.searchfactor;

nvars = length(x);

sigma = optimState.C';
sigma = MeshSize*SearchFactor*sigma;    % Rescale by current scale

N = options.Nsearch;
xs = bsxfun(@plus, x, randn(N,nvars)*sigma);

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
% IN ORDER TO SIMULATE THE PROGRAM:
%   A) FIRST, LOAD A ROBOT
%       robot = load_robot('abb','irb140');
%   B) NEXT, LOAD SOME EQUIPMENT.
%       robot.equipment = load_robot('equipment','tables/table_small');
%       OR
%       robot.equipment = load_robot('equipment','bumper_cutting');
%   C) NOW, LOAD AN END TOOL
%       robot.tool= load_robot('equipment','end_tools/parallel_gripper_0');
%   D) FINALLY, LOAD A PIECE TO GRAB BY THE ROBOT
%       robot.piece=load_robot('equipment','cylinders/cylinder_tiny');
%
%   E) IF NECESSARY, CHANGE THE POSITION AND ORIENTATION OF THE ROBOT'S
%   BASE
%       robot.piece.T0= [1 0 0 -0.35;
%                        0 1 0 -0.55;
%                        0 0 1 0.2;
%                        0 0 0 1];
%
% during the simulation, call simulation_open_tool; to open the tool and
% simulation_close_tool; to close it.
% To grip the piece, call simulation_grip_piece; and
% simulation_release_piece to release it.
% The call to each function must be correct, thus, typically the correct
% sequence is:
% simulation_open_tool;
% approach the piece to grab.
% simulation_close_tool;
% simulation_grip_piece; --> the piece will be drawn with the robot
% move to a different place
% simulation_open_tool;
% simulation_release_piece

function test_1

global RT_tp1 RT_tp2 TD_tool0

RT_tp1=[[0.4000, 0.3000, 0.7000],[1.0000, -0.0000, -0.0000, -0.0000], [0, -1, -1, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
RT_tp2=[[0.1000, -0.4000, 0.7000],[0.0000, -0.0000, 0.0000, 1.0000], [-1, -1, -2, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];

%definition of the end effector.
TD_tool0=[1,[[0,0,0],[1,0,0,0]],[0,[0,0,0],[1,0,0,0],0,0,0]];

%finally, call the main function
main
end

function main
global RT_tp1 RT_tp2 TD_tool0

MoveJ(RT_tp1, 'vmax' , 'fine' , TD_tool0, 'wobj0');
MoveJ(RT_tp2, 'vmax' , 'fine' , TD_tool0, 'wobj0');
MoveL(RT_tp1, 'vmax' , 'fine' , TD_tool0, 'wobj0');
MoveL(RT_tp2, 'vmax' , 'fine' , TD_tool0, 'wobj0');

end

repo: ARTE


text:
%  Figure 3.2      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%script to generate Fig. 3.2
%% fig3_02.m    Example 3.4
clf;
k=1;
num=1;                        % form numerator
den=[1 k];                    % form denominator
% sinusoidal input signal
deltaT = 0.001;
t=0:deltaT:10;                % form time vector
u=sin(10*(t));                % form input
sys=tf(num,den);              % form system
[y]=lsim(sys,u,t);            % linear simulation
%plot response
plot(t,y);
xlabel('Time (sec)');
ylabel('Output');
title('Fig. 3.2 (a): transient response');
grid;
pause;
hold on;
y1=(10/101)*exp(-t);
phi=atan(-10);
y2=(1/sqrt(101))*sin(10*t+phi);
plot(t,y1,t,y2,t,y1+y2);
pause;
hold off;
ii=[9001:10001];
plot(t(ii),y(ii),t(ii),u(ii));
xlabel('Time (sec)');
ylabel('Output, input');
title('Fig. 3.2 (b): Steady-state response');
grid;
text(9.4,0.65,'u(t)');
text(9.24,0.12,'y(t)');


text:
function y = daxpy ( n, sa, x, incx, y, incy )

%*****************************************************************************80
%
%% DAXPY adds a constant times one vector to another.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real SA, the multiplier.
%
%    Input, real X(*), the vector to be scaled and added to Y.
%
%    Input, integer INCX, the increment between successive entries of X.
%
%    Input/output, real Y(*), the vector to which a
%    multiple of X is to be added.
%
%    Input, integer INCY, the increment between successive entries of Y.
%
  y(1:incy:1+(n-1)*incy) = y(1:incy:1+(n-1)*incy) + sa * x(1:incx:1+(n-1)*incx);

  return
end


text:
function test_bug1646

% MEM 6gb
% WALLTIME 00:15:00
% DEPENDENCY ft_prepare_mesh ft_datatype_segmentation

% the purpose of this test script is to ensure that the new implementation
% of ft_prepare_mesh, which is a merger between the old ft_prepare_mesh and
% ft_prepare_mesh_new, works as expected on all normal cases. What I have
% done is to search for as many as possible other tets cases and replicated
% them here.

% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=105 is obsolete
% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=115 is obsolete
% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=964 is obsolete

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

% it is not supposed to work with a 3-D array as input
try
  cfg = [];
  bnd = ft_prepare_mesh(cfg, randn(30,30,30));
  status = false;
catch
  status = true;
end
if ~status
  error('invalid input was not dealt with')
end

% it is supposed to wotk with a headmodel as input
headmodel = [];
headmodel.r = 12;
headmodel.o = [0 0 4];
headmodel.unit = 'cm';

cfg = [];
cfg.numvertices = 100;
bnd = ft_prepare_mesh(cfg, headmodel);

cfg = [];
bnd = ft_prepare_mesh(cfg, headmodel);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the following piece of test code relates loosely to
% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=222
% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=373
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cfg = [];
cfg.headshape = dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.shape');
bnd = ft_prepare_mesh(cfg);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is the same test data as used for bug 1652
% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1652
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

datadir = fullfile(dccnpath('/home/common/matlab/fieldtrip/data/test/bug1646'));

load(fullfile(datadir,'seg1'));
load(fullfile(datadir,'seg2'));
load(fullfile(datadir,'seg3'));
load(fullfile(datadir,'seg4'));
load(fullfile(datadir,'seg5'));
load(fullfile(datadir,'seg6'));

atlas = ft_read_atlas(dccnpath('/home/common/matlab/fieldtrip/template/atlas/afni/TTatlas+tlrc.BRIK'));
mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer/Subject01.mri'));

assert(ft_datatype(seg1, 'segmentation'));
assert(ft_datatype(seg2, 'segmentation'));
assert(ft_datatype(seg3, 'segmentation'));
assert(ft_datatype(seg4, 'segmentation'));
assert(ft_datatype(seg5, 'segmentation'));
assert(ft_datatype(seg6, 'segmentation'));
assert(ft_datatype(atlas, 'segmentation'));
assert(ft_datatype(mri, 'volume'));

cfg = [];
cfg.numvertices = 642;
cfg.tissue = 'brain';
bnd = ft_prepare_mesh(cfg, seg1);
assert(isequal(bnd(1).unit, seg1.unit));

% this is a weird case, because seg2 is not indexed but tpm
% the output depends on the order of the cfg.tissue
cfg = [];
cfg.numvertices = [1000, 100, 10];
cfg.tissue = [1 2 3];
bndA = ft_prepare_mesh(cfg, seg2);
cfg.numvertices = [10, 100, 1000];
cfg.tissue = [3 2 1];
bndB = ft_prepare_mesh(cfg, seg2);
assert(isequal(bndA(1), bndB(3)));

cfg = [];
cfg.numvertices = [1000, 100, 10];
cfg.tissue = {'brain', 'skull', 'scalp'};
bndA = ft_prepare_mesh(cfg, seg2);
assert(isequal(bndA(1).unit, seg2.unit));
assert(size(bndA(1).pos,1)==1000);
assert(size(bndA(2).pos,1)==100);
assert(size(bndA(3).pos,1)==10);
% check that the output has the desired order (as specified in the cfg)
cfg.numvertices = [10, 100, 1000];
cfg.tissue = {'scalp', 'skull', 'brain'};
bndB = ft_prepare_mesh(cfg, seg2);
assert(isequal(bndA(1), bndB(3)));

cfg = [];
cfg.numvertices = [1000, 100, 10];
cfg.tissue = {'brain', 'skull', 'scalp'};
bnd = ft_prepare_mesh(cfg, seg3);
assert(isequal(bnd(1).unit, seg3.unit));
assert(size(bnd(1).pos,1)==1000);
assert(size(bnd(2).pos,1)==100);
assert(size(bnd(3).pos,1)==10);

cfg = [];
cfg.numvertices = [1000, 100, 10];
cfg.tissue = {'brain', 'skull', 'scalp'};
bnd = ft_prepare_mesh(cfg, seg4);
assert(isequal(bnd(1).unit, seg4.unit));
assert(size(bnd(1).pos,1)==1000);
assert(size(bnd(2).pos,1)==100);
assert(size(bnd(3).pos,1)==10);

cfg = [];
cfg.numvertices = [1000, 100, 10];
cfg.tissue = [3 2 1];
bnd = ft_prepare_mesh(cfg, seg5);
assert(isequal(bnd(1).unit, seg5.unit));
assert(size(bnd(1).pos,1)==1000);
assert(size(bnd(2).pos,1)==100);
assert(size(bnd(3).pos,1)==10);

cfg = [];
cfg.numvertices = 642;
cfg.tissue = 'brain';
bnd = ft_prepare_mesh(cfg, seg6);
assert(isequal(bnd(1).unit, seg6.unit));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% addressing further issues raised in bug1646, on cfg.tissue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% If the input is an indexed representation and  cfg.tissue is equal to
% setdiff(unique(mri.seg(:)),0) then there is no problem.

% seg5 is indexed segmentation
cfg = [];
cfg.tissue = setdiff(unique(seg5.seg(:)),0);
cfg.numvertices = 1000;
bnd = ft_prepare_mesh(cfg, seg5);

% "If the input is an indexed representation and cfg.tissue is empty, then
% cfg.tissue should be set to unique(mri.seg(:))"

cfg = [];
cfg.numvertices = [1000 1000 1000];
cfg.tissue = [];
bnd = ft_prepare_mesh(cfg, seg5);

cfg = [];
cfg.numvertices = [1000 1000 1000];
bnd = ft_prepare_mesh(cfg, seg5);

cfg = [];
cfg.numvertices = [1000]; % should still be converted to vector Nx3 since seg5 has 3 indices
cfg.tissue = [];
bnd = ft_prepare_mesh(cfg, seg5);

% "If the input is a tpm representation then cfg.tissue should be a string
% pointing to the field."
% Both seg1 and seg2 are tpm, but only seg2 has 'brain'.

cfg = [];
cfg.numvertices = [800];
cfg.tissue = 'brain'; %this is exception in that .brain DNE in seg1
bnd = ft_prepare_mesh(cfg, seg1);

try
  cfg = [];
  cfg.numvertices = [800];
  cfg.tissue = 'randomfieldname';
  bnd = ft_prepare_mesh(cfg, seg1);
  success = true;
catch
  success = false;
end
if success
  error('randomfieldname should not work for tpm type')
end

try
  cfg = [];
  cfg.numvertices = [800];
  cfg.tissue = 'csf'; %this is exception in that .brain DNE in seg1
  bnd = ft_prepare_mesh(cfg, seg2);
  success = true;
catch
  success = false;
end
if success
  error('randomfieldname should not work for tpm type')
end

cfg = [];
cfg.numvertices = [800];
cfg.tissue = 1;
bnd = ft_prepare_mesh(cfg, seg1);

cfg = [];
cfg.numvertices = [800];
cfg.tissue = 1;
bnd = ft_prepare_mesh(cfg, seg2);

try
  cfg = [];
  cfg.numvertices = [800];
  cfg.tissue = 4;
  bnd = ft_prepare_mesh(cfg, seg2);
  success = true;
catch
  success = false;
end
if success
  error('too large of index for cfg.tissue should not work for tpm type')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this shares the test data with bug 1651
% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1651
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug1651.mat'));

cfg = [];
cfg.tissue = {'brain', 'skull', 'scalp'};
cfg.numvertices = [1000 2000 3000];
bnd = ft_prepare_mesh(cfg, seg2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=997
% this one does not state an explicit way of reproducing the problems, it rather refers to 937
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=937
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

mri = [];
mri.anatomy = bkgrnd;
mri.sphere1 = bkgrnd==1;
mri.sphere2 = bkgrnd==2;
mri.sphere3 = bkgrnd==3;
mri.transform = eye(4);
mri.dim = size(bkgrnd);

mri.anatomy = mri.anatomy>0;

cfg = [];
cfg.tissue = {'sphere1' 'sphere2' 'sphere3'};
cfg.numvertices = 1000;
bnd = ft_prepare_mesh(cfg, mri);


text:
function inside = checkInside(oR,q,varargin)
% check for points to be inside the orientation region

% avoid q beeing orientation
q = quaternion(q);

if isempty(oR.N), inside = true(size(q)); return; end
if isempty(q), inside = false(size(q)); return; end

% verify all conditions are satisfies
inside1 = dot_outer(oR.N,q)<1e-3;
inside2 = dot_outer(oR.N,q)>-1e-3;

% either q or -q needs to satisfy the condition
inside = reshape(all(inside1,1) | all(inside2,1),size(q));

end

author: mtex-toolbox

repo: mtex

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


text:
% VL_LBP  Local Binary Patterns
%   F = VL_LBP(IM, CELLSIZE) computes the Local Binary Pattern (LBP)
%   features for image I.
%
%   IM is divided in cells of size CELLSIZE. F is a three-dimensional
%   array containing one histograms of quantized LBP features per
%   cell. The witdh of F is FLOOR(WIDTH/CELLSIZE), where WIDTH is the
%   width of the image. The same for the height. The third dimension
%   is 58.
%
%   See also: VL_HELP().

% AUTORIGHTS