text:
function varargout = gsynsq(varargin)
% function hgui = gsynsq()
%
% Runs the Synchrosqueezing toolbox GUI.
% Does not require input parameters.  Returns a handle to the GUI.
%
% This GUI tool can be used for speedy/rough Synchrosqueezing
% analysis, filtering, component extraction of signals.
%
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/)
%---------------------------------------------------------------------------------

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

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

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

% Choose default command line output for gsynsq
hs.output = hObj;

% This sets up the initial plot - only do when we are invisible
% so window can get raised using gsynsq.
if strcmp(get(hObj,'Visible'),'off')
    set(hs.axesx, 'Visible', 'off');
    set(hs.axesTx, 'Visible', 'off');
    set(hs.axesrx, 'Visible', 'off');

    % Set up zoom
    zh = zoom;

    setAllowAxesZoom(zh, hs.axesx, true);
    setAllowAxesZoom(zh, hs.axesTx, true);
    setAllowAxesZoom(zh, hs.axesrx, true);

    % Link the axes w.r.t. zoom (just the t-axis component)
    linkaxes([hs.axesx hs.axesTx hs.axesrx], 'x');

    % Menu bar with "reset axes" does not work with linked axes
    set(zh, 'RightClickAction', 'InverseZoom');

    % Append any past ActionPostCallBack to run first
    %set(zh, 'ActionPostCallback', @postZoomUpdate);

    % Allow zoom
    set(zh, 'Enable', 'on');

    hs.zh = zh;

    % Initialize internal data structure placeholders
    hs = clear_data(hs);
end

% Update hs structure
guidata(hObj, hs);

% UIWAIT makes gsynsq wait for user response (see UIRESUME)
% uiwait(hs.figure1);

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

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

% --------------------------------------------------------------------
function FileMenu_Callback(hObj, eventdata, hs)
% hObj    handle to FileMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)

% --------------------------------------------------------------------
function OpenMenuItem_Callback(hObj, eventdata, hs)
% hObj    handle to OpenMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
file = uigetfile('*.mat');
if ~isequal(file, 0)
    try
        tx = load(file);
        assert(isfield(tx, 't'));
        assert(isfield(tx, 'x'));
        assert(length(tx.t) == length(tx.x));
    catch ME
        errordlg(sprintf('File load error: %s', ME.message), ...
                 'File Load Error');
        return;
    end

    load_data(hObj, eventdata, hs, file, tx);
end

function load_data(hObj, eventdata, hs, file, tx)
% Do something here about old data, if any?
if (hs.active && ~hs.saved)
end

% Load up the data.  First clear old data
hs = clear_data(hs);
hs.data.t = tx.t;
hs.data.x = tx.x - mean(tx.x);

[fpath, fname] = fileparts(file);
hs.prefix = fullfile(fpath, fname);

% Plot axes
plot(hs.axesx, hs.data.t, hs.data.x);
set(hs.axesx, 'Visible', 'on');
grid(hs.axesx, 'on');
axis(hs.axesx,'tight');
setAxesZoomMotion(hs.zh, hs.axesx, 'horizontal');
set(hs.axesx,'FontSize',10);

% Hide other axes and delete colorbars, if any
axes(hs.axesTx);
colorbar('delete');

cla(hs.axesTx);
set(hs.axesTx, 'Visible', 'off');
cla(hs.axesrx);
set(hs.axesrx, 'Visible', 'off');

% Save app data
guidata(hObj, hs);

% --------------------------------------------------------------------
function PrintMenuItem_Callback(hObj, eventdata, hs)
% hObj    handle to PrintMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
printdlg(hs.figure1)

% --------------------------------------------------------------------
function ExitMenuItem_Callback(hObj, eventdata, hs)
% hObj    handle to CloseMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
figure1_CloseRequestFcn(hObj, eventdata, hs);

% --- Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObj, eventdata, hs)
% hObj    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)

% If the analysis has been modified but not been saved
if (hs.active && ~hs.saved)
    selection = questdlg(['Unsaved Analysis.  Close ' get(hs.figure1,'Name') '?'],...
        ['Close ' get(hs.figure1,'Name') '...'],...
        'Yes','No','Yes');
    if strcmp(selection,'No')
        return;
    end
end

% Closes the figure
delete(hs.output);

% Initializes / clears data from last job from hs
function hs = clear_data(hs)
% Is there an analysis active?
hs.active = 0;
% Has it been saved since the last change?
hs.saved = 0;

% What is the file prefix of the current analysis?
hs.prefix = '';

% Default parameters for the transforms
hs.padtype = 'symmetric';
% Number of voices
hs.nv = 32;
% Wavelet type and other options
hs.opt = struct('disp', 1, 'bd', 1, 'gamma', sqrt(eps), 'type', 'morlet');

% Placeholders for the input data
hs.data = [];

% Placeholders for the results of the analysis part
hs.data.as = [];
hs.data.fs = [];
hs.data.Wx = [];
hs.data.Tx = [];

% Placeholders for the results of the synthesis part
hs.data.rx = [];

% Placeholder for the results of synthesis after removing padding
%hs.data.rx0 = [];

% Placeholder for the extracted curves
hs.data.Cs = [];
hs.data.xrs = [];

% Attenuation filters
if isfield(hs, 'filters')
    for hi = 1:length(hs.filters)
        delete(hs.filters(hi).f);
    end
end
hs.filters = [];

% --------------------------------------------------------------------
function openpush_ClickedCallback(hObj, eventdata, hs)
% hObj    handle to openpush (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
OpenMenuItem_Callback(hObj, eventdata, hs);

% --------------------------------------------------------------------
function savepush_ClickedCallback(hObj, eventdata, hs)
% hObj    handle to savepush (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
SaveMenuItem_Callback(hObj, eventdata, hs);

% --------------------------------------------------------------------
function zoomtoggle_ClickedCallback(hObj, eventdata, hs)
% hObj    handle to zoomtoggle (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)

% ZOOM

% --------------------------------------------------------------------
function zoomtoggle_OffCallback(hObj, eventdata, hs)
% hObj    handle to zoomtoggle (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)

% Disable zoom
zoom_disable(hs);

function zoom_disable(hs)
set(hs.zoomtoggle, 'State', 'off');
set(hs.ZoomMenuItem, 'Checked', 'off');
set(hs.zh, 'Enable', 'off');

% --------------------------------------------------------------------
function zoomtoggle_OnCallback(hObj, eventdata, hs)
% hObj    handle to zoomtoggle (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)

% Allow zoom
%if strcmp(get(hObj,'Visible'),'off')
zoom_enable(hs);
%end

function zoom_enable(hs)
set(hs.zoomtoggle, 'State', 'on');
set(hs.zh, 'Enable', 'on');
set(hs.ZoomMenuItem, 'Checked', 'on');

%%% NOT USED %%%
function postZoomUpdate(hObj, eventdata)
% Get the new horizontal values of the updated axes
XLim = get(eventdata.Axes, 'XLim');

% Update the values for each of the axes
hs = guidata(hObj);

XLim(1) = max(XLim(1), min(hs.data.t));
XLim(2) = min(XLim(2), max(hs.data.t));
set(eventdata.Axes, 'XLim', XLim);

if (eventdata.Axes ~= hs.axesTx), set(hs.axesTx, 'XLim', XLim); end
if (eventdata.Axes ~= hs.axesx), set(hs.axesx, 'XLim', XLim); end
if (eventdata.Axes ~= hs.axesrx), set(hs.axesrx, 'XLim', XLim); end

% --------------------------------------------------------------------
function SaveMenuItem_Callback(hObj, eventdata, hs)
% hObj    handle to SaveMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
if (~hs.active || hs.saved)
    return;
end

file = uiputfile(sprintf('%s_r.mat', hs.prefix));
if ~isequal(file, 0)
    try
        data = hs.data;
        save(file, '-struct', 'data');
    catch ME
        errordlg(sprintf('File save error: %s', ME.message), ...
            'File Save Error');
        return;
    end

    hs.saved = 1;
    guidata(hObj, hs);
end

% --------------------------------------------------------------------
function ExportMenuItem_Callback(hObj, eventdata, hs)
% hObj    handle to ExportMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
if (~hs.active), return; end
labels = {'Export t?', 'Export x(t)?', 'Export as?', 'Export Wx(a,t)?', ...
          'Export fs?', 'Export Tx(f,t)?', 'Export rx(t) (reconstructed)?', ...
          'Export filters?', 'Export identified contours?', ...
          'Export extracted components?'};
varnames = {'t', 'x', 'as', 'Wx', 'fs', 'Tx', 'rx', ...
            'filters', 'Cs', 'xrs'};
filters = hs.filters;
for fi=1:length(filters)
    filters(fi).bbox = filters(fi).f.getPosition();
    filters(fi).mask = filters(fi).f.createMask();
    % bbox = [tmin tmax fmin fmax]
    % filters(fi).bbox = [filters(fi).bbox([1 3])];
    tmin = filters(fi).bbox(1);
    tmax = tmin + filters(fi).bbox(3);
    flogmin = filters(fi).bbox(2);
    flogmax = flogmin + filters(fi).bbox(4);
    filters(fi).bbox = [tmin tmax 2^flogmin 2^flogmax];
end
if length(filters)>0, filters = rmfield(filters, 'f'); end
vars = {hs.data.t, hs.data.x, hs.data.as, ...
        hs.data.Wx, hs.data.fs, hs.data.Tx, hs.data.rx, ...
        filters, hs.data.Cs, hs.data.xrs};
selected = logical([0 0 0 0 0 0 0 0 0 0]);
export2wsdlg(labels, varnames, vars, 'Export to Workspace', selected);

% TODO - export to a variable with proper prefix and name of choice

% --------------------------------------------------------------------
function AnalysisMenu_Callback(hObj, eventdata, hs)
% hObj    handle to AnalysisMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)

% --------------------------------------------------------------------
function AnalysisMenuItem_Callback(hObj, eventdata, hs)
% hObj    handle to AnalysisMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
run_analysis(hObj, eventdata, hs);

function run_analysis(hObj, eventdata, hs)
if (isempty(hs.data))
    errordlg('Load a data file first', 'Error');
    return;
end

prompt = {'Wavelet [morlet,bump,hhhat,...]', ...
          'Number of Voices (nv)', ...
          'Pad type [symmetric, circular, replicate]', ...
          'Gamma (\gamma)'};
name = 'Synchrosqueezing Analysis Parameters';
answers = {hs.opt.type, num2str(hs.nv), hs.padtype, num2str(hs.opt.gamma)};
idata = inputdlg(prompt, name, 1, answers, struct('WindowStyle', 'normal'));
if isempty(idata)
    return
end

hs.opt.type = idata{1};
hs.nv = str2num(idata{2});
hs.padtype = idata{3};
hs.opt.gamma = str2num(idata{4});

[hs.data.Tx, hs.data.fs, hs.data.Wx, hs.data.as] = ...
    synsq_cwt_fw(hs.data.t, hs.data.x, hs.nv, hs.opt);

axes(hs.axesTx);
cla(hs.axesTx);
tplot(hs.data.Tx, hs.data.t, hs.data.fs, hs.opt);
set(hs.axesTx, 'Visible', 'on');
grid(hs.axesTx, 'on');
setAxesZoomMotion(hs.zh, hs.axesTx, 'vertical');
set(hs.axesTx,'FontSize',10);
hs.filters = [];

hs.active = 1;

% Store application parameters
guidata(hObj, hs);

% --------------------------------------------------------------------
function ReconMenuItem_Callback(hObj, eventdata, hs)
% hObj    handle to ReconMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
run_recon(hObj, eventdata, hs);

function run_recon(hObj, eventdata, hs);

if (~hs.active)
    errordlg('Run Analysis First', 'Error');
    return;
end

Tx = filter_Tx(hs);

hs.data.rx = synsq_cwt_iw(Tx, hs.data.fs, hs.opt);

plot(hs.axesrx, hs.data.t, hs.data.rx);
grid(hs.axesrx,'on');
set(hs.axesrx, 'Visible', 'on');
set(hs.axesrx,'FontSize',10);
axis(hs.axesrx, [min(hs.data.t) max(hs.data.t) min(hs.data.x) max(hs.data.x)]);

XLim = get(hs.axesx, 'XLim'); set(hs.axesrx, 'XLim', XLim);
setAxesZoomMotion(hs.zh, hs.axesrx, 'horizontal');
%zoom_enable(hs);

guidata(hObj, hs);

function Tx = filter_Tx(hs)
% Logic to see if any attenuation and/or passband filters are active
Tx = hs.data.Tx;
if ~isempty(hs.filters)
    % Step 1: combine all pass-band filters
    maskp = zeros(size(hs.data.Tx));
    maska = ones(size(hs.data.Tx));
    nmp = 0;
    for fi=1:length(hs.filters)
        mask = hs.filters(fi).f.createMask();
        if strcmpi(hs.filters(fi).type, 'atten')
            maska = maska & ~mask;
        else
            maskp = maskp | mask;
            nmp = nmp + 1;
        end
    end

    % Mask via the filters
    if (nmp > 0)
        Tx = Tx .* maskp .* maska;
    else
        Tx = Tx .* maska;
    end
end

% --------------------------------------------------------------------
function FiltersMenu_Callback(hObj, eventdata, hs)
% hObj    handle to FiltersMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)

% --------------------------------------------------------------------
function RectAtten_Callback(hObj, eventdata, hs)
% hObj    handle to RectAtten (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
filter_enable(hObj, hs, 'recta')

% --------------------------------------------------------------------
function filter_enable(hObj, hs, type)

% Check to make sure analysis is active
if (~hs.active)
    return;
end

% Disable zoom icon
set(hs.zoomtoggle, 'State', 'off');

% Start filter selection - this disabled zoom
fhs = filter_construct(hs, type);

% Make sure that it's within the proper constraints at first
for fhi=1:length(fhs)
    fh = fhs(fhi);
    fh.setConstrainedPosition(fh.getPosition());

    if strcmpi(type(end), 'a')
        hs.filters = [hs.filters struct('type', 'atten', 'f', fh)];
    else
        hs.filters = [hs.filters struct('type', 'pass', 'f', fh)];
    end
end

% Save
guidata(hObj, hs);

function h = filter_construct(hs, type)
switch type
  case 'recta',
    h = imrect(hs.axesTx);
  case 'ellipsea',
    h = imellipse(hs.axesTx);
  case 'polya',
    h = impoly(hs.axesTx);
  case 'rectp',
    h = imrect(hs.axesTx);
  case 'bda',
    t = hs.data.t(:);
    fs = [hs.data.fs(1), hs.data.fs(2:hs.nv:end-1), hs.data.fs(end)]; % Subsample
    [Lbdi, Rbdi] = synsq_bd_loc(t, fs);
    fsbd = [fs(1); fs(:); fs(end)];
    Lbdi = [1; Lbdi(:); 1];
    Rbdi = [length(t); Rbdi(:); length(t)];

    hL = impoly(hs.axesTx, [t(Lbdi) log2(fsbd)]);
    hR = impoly(hs.axesTx, [t(Rbdi) log2(fsbd)]);
    h = [hL hR];
  otherwise,
    error('Unknown filter type');
end

for j=1:length(h)
    hj = h(j);

    cstrstr = class(hj);

    % Set constraints
    constrfcn = makeConstrainToRectFcn(cstrstr, ...
      [min(hs.data.t)-sqrt(eps) max(hs.data.t)], ...
      [min(log2(hs.data.fs))-sqrt(eps) max(log2(hs.data.fs))]);
    hj.setPositionConstraintFcn(constrfcn);

    % Set color
    if strcmpi(type(end), 'a')
        hj.setColor('r');
    else
        hj.setColor('g');
    end
end

% --- Executes on mouse press over axes background.
function axesTx_ButtonDownFcn(hObj, eventdata, hs)
% hObj    handle to axesTx (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)

% Is not called due to axis ownership by zoom object

% --------------------------------------------------------------------
function EllipseAtten_Callback(hObj, eventdata, hs)
% hObj    handle to EllipseAtten (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
filter_enable(hObj, hs, 'ellipsea')

% --------------------------------------------------------------------
function PolyAtten_Callback(hObj, eventdata, hs)
% hObj    handle to PolyAtten (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
filter_enable(hObj, hs, 'polya')

% --------------------------------------------------------------------
function filterselector_ClickedCallback(hObj, eventdata, hs)
% hObj    handle to filterselector (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)

% FILTER SELECTOR

% --------------------------------------------------------------------
function filterselector_OffCallback(hObj, eventdata, hs)
% hObj    handle to filterselector (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)

% FILTER SELECTOR

% --------------------------------------------------------------------
function filterselector_OnCallback(hObj, eventdata, hs)
% hObj    handle to filterselector (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)

% FILTER SELECTOR

% --------------------------------------------------------------------
function rectattenpush_ClickedCallback(hObj, eventdata, hs)
% hObj    handle to rectattenpush (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
filter_enable(hObj, hs, 'recta')

% --------------------------------------------------------------------
function polyattenpush_ClickedCallback(hObj, eventdata, hs)
% hObj    handle to polyattenpush (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
filter_enable(hObj, hs, 'polya')

% --------------------------------------------------------------------
function ellipseattenpush_ClickedCallback(hObj, eventdata, hs)
% hObj    handle to ellipseattenpush (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
filter_enable(hObj, hs, 'ellipsea')

% --------------------------------------------------------------------
function analysispush_ClickedCallback(hObj, eventdata, hs)
% hObj    handle to analysispush (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
run_analysis(hObj, eventdata, hs);

% --------------------------------------------------------------------
function reconpush_ClickedCallback(hObj, eventdata, hs)
% hObj    handle to reconpush (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% hs    structure with hs and user data (see GUIDATA)
run_recon(hObj, eventdata, hs);

% --------------------------------------------------------------------
function ZoomMenuItem_Callback(hObject, eventdata, handles)
% hObject    handle to ZoomMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
if strcmpi(get(handles.zh, 'Enable'), 'on'),
    zoomtoggle_OffCallback(hObject, eventdata, handles);
else
    zoomtoggle_OnCallback(hObject, eventdata, handles);
end

% --------------------------------------------------------------------
function ToolsMenu_Callback(hObject, eventdata, handles)
% hObject    handle to ToolsMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function ImportMenuItem_Callback(hObject, eventdata, handles)
% hObject    handle to ImportMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

prompts = {'Import t variable: ', 'Import x(t) variable: ', 'Prefix: '};
name = 'Import signal from workspace';
danswers = {'t', 'x', 'ws'};
idata = inputdlg(prompts, name, 1, danswers, struct('WindowStyle', 'normal'));

% Cancelled
if isempty(idata),
    return;
end

try
    tx = struct('t', evalin('base', idata{1}), ...
        'x', evalin('base', idata{2}));
    assert(length(tx.t) == length(tx.x));
catch ME
    errordlg(sprintf('Import error: %s', ME.message), ...
        'Import Error');
        return;
end

% Set some file prefix in the current directory
file = fullfile(pwd, idata{3});

load_data(hObject, eventdata, handles, file, tx);

% --------------------------------------------------------------------
function rectpasspush_ClickedCallback(hObj, eventdata, hs)
% hObject    handle to rectpasspush (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
filter_enable(hObj, hs, 'rectp');

% --------------------------------------------------------------------
function RectPass_Callback(hObj, eventdata, hs)
% hObject    handle to RectPass (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
filter_enable(hObj, hs, 'rectp');

% --------------------------------------------------------------------
function ExtractMenuItem_Callback(hObject, eventdata, handles)
% hObject    handle to ExtractMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
extract_contour(hObject, handles);

% Extract a contour
function extract_contour(hObj, hs)

if (~hs.active)
    errordlg('Run Analysis First', 'Error');
    return;
end

prompt = {'Number of Curves to Extract', ...
          'Lambda (smoothness penalty)', ...
         'Extraction window (n_w)'};
name = 'Contour Extraction Parameters';
numlines = 1;
answers = {'2', '1e4', num2str(hs.nv/2)};

idata = inputdlg(prompt, name, numlines, answers, struct('WindowStyle', 'normal'));

% User cancelled
if isempty(idata),
    return;
end

nc = str2num(idata{1});
lambda = str2num(idata{2});
nw = str2num(idata{3});

Tx = filter_Tx(hs);
[na, N] = size(Tx);

% Extract curves for plotting
[hs.data.Cs,Es] = curve_ext_multi(Tx, log2(hs.data.fs), nc, lambda, nw);

% Extract functions for export
hs.data.xrs = curve_ext_recon(Tx, hs.data.fs, hs.data.Cs, hs.opt, nw);

% Save
guidata(hObj, hs);

figure;
hc = plot_ext_curves(hs.data.t, hs.data.x, ...
                     Tx, hs.data.fs, hs.data.Cs, Es, hs.opt);

% --------------------------------------------------------------------
function extractpush_ClickedCallback(hObject, eventdata, handles)
% hObject    handle to extractpush (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
extract_contour(hObject, handles);

% --------------------------------------------------------------------
function AttenBd_Callback(hObj, eventdata, hs)
% hObject    handle to AttenBd (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
filter_enable(hObj, hs, 'bda');

% --------------------------------------------------------------------
function attenbdpush_ClickedCallback(hObj, eventdata, hs)
% hObject    handle to attenbdpush (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
filter_enable(hObj, hs, 'bda');

author: HeLiangHIT


text:
function clkeysel()
% provides a interface for some selection options by a polygon
% works on the equivalent events for a cluster catalog(complete,swarms,etc)
% selects clusters which equivalent events are inside selection area
%
%A.Allmann
 % turned into function by Celso G Reyes 2017

ZG=ZmapGlobal.Data; % used by get_zmap_globals
global x y clu mess plot1_h plot2_h clust file1
global  h5 xcordinate ycordinate newclcat clsel
global equi_button backbgevent original backcat backequi
global decc

xcordinate=0;
ycordinate=0;
axes(h5)
x = [];
y = [];

%n = 0;

figure(mess);
set(gcf,'visible','off')
clf;
cla;
set(gca,'visible','off');
set(gcf, 'Name','Polygon Input Parameters');

%creates dialog box to input some parameters
%

inp1_field=uicontrol('Style','edit',...
    'Position',[.70 .60 .17 .10],...
    'Units','normalized','String',num2str(xcordinate),...
    'callback',@callbackfun_001);

inp2_field=uicontrol('Style','edit',...
    'Position',[.70 .40 .17 .10],...
    'Units','normalized','String',num2str(ycordinate),...
    'callback',@callbackfun_002);

more_button=uicontrol('Style','Pushbutton',...
    'Position', [.60 .05 .15 .15],...
    'Units','normalized',...
    'callback',@callbackfun_003,...
    'String','More');
last_button=uicontrol('Style','Pushbutton',...
    'Position',[.40 .05 .15 .15],...
    'Units','normalized',...
    'Callback',@(~,~) clpickp('LAST'),...
    'String','Last');

mouse_button=uicontrol('Style','Pushbutton',...
    'Position',[.20 .05 .15 .15],...
    'Units','normalized',...
    'Callback',@(~,~)clpickp('MOUSE'),...
    'String','Mouse');

load_button=uicontrol('Style','Pushbutton',...
    'Position',[.80 .05 .15 .15],...
    'Units','normalized',...
    'Callback',@(~,~)clpickp('LOAD'),...
    'String','load');
cancel_button=uicontrol('Style','Pushbutton',...
    'Position',[.05 .80 .15 .15],...
    'Units','normalized',...
    'callback',@callbackfun_004,...
    'String','Cancel');
txt1 = text(...
    'Position',[0. 0.65 0 ],...
    'FontSize',12,...
    'String','Longitude:');
txt2 = text(...
    'Position',[0. 0.45 0 ],...
    'FontSize',12,...
    'String','Latitude:');

set(mess,'visible','on')

function callbackfun_001(mysrc,myevt)

  callback_tracker(mysrc,myevt,mfilename('fullpath'));
  xcordinate=str2double(inp1_field.String);
  inp1_field.String=num2str(xcordinate);
end

function callbackfun_002(mysrc,myevt)

  callback_tracker(mysrc,myevt,mfilename('fullpath'));
  ycordinate=str2double(inp2_field.String);
  inp2_field.String=num2str(ycordinate);
end

function callbackfun_003(mysrc,myevt)

  callback_tracker(mysrc,myevt,mfilename('fullpath'));
  set(mouse_button,'visible','off');
  set(load_button,'visible','off');
  clpickp('MORE');
end

function callbackfun_004(mysrc,myevt)

  callback_tracker(mysrc,myevt,mfilename('fullpath'));
  close
end

end

author: CelsoReyes


text:
function hObj = findobj2(h, propname, propvalue, options)

% Find graphic object handle given the parent (or grandparent, etc) handle
% and property name/value pair

if ~exist('options','var')
    options = '';
end

hObj = [];
hc = get(h, 'children');
for ii = 1:length(hc)
    try
        if eval( sprintf('strcmp(get(hc(ii), ''%s''), ''%s'')', propname, propvalue) )
            hObj = hc(ii);
            return;
        elseif ~optionExist(options, 'flat')
            hObj = findobj2(hc(ii), propname, propvalue);
            if ~isempty(hObj)
                return;
            end
        end
    catch
        continue
    end
end

author: BUNPC


text:
function view = motionCompNestaresWithin1st(view, tgtScans, baseScan, baseFrame, smoothing);
%
% view = motionCompNestaresFull(view, <tgtScans=all>, <baseScan=1>, , <baseFrame=1>, <temporalSmoothing = 1 frame>);
%
% Call both between and within scan motion compensation, back-to-back, on
% an inplane view. Right now just a simple wrapper.
%
% Note that these are the older, rigid-body only motion compensation tools
% based around the code originally developed by Nestares and Heeger.
%
%
% DY 09/2008 based on motionCompNestaresFull
if notDefined('view'), view = getSelectedInplane; end
if notDefined('tgtScans'), tgtScans = 1:numScans(view); end
if notDefined('baseScan'), baseScan = 1; end
if notDefined('baseFrame'), baseFrame = 1; end
if notDefined('smoothing'), smoothing = 1; end

if ~isequal(view.viewType, 'Inplane')
    myErrorDlg('Can only run motion compensation on Inplane data.');
end

% first do within scans motion compensation:
newDataType = 'MotionComp';
view = motionCompSelScan(view, newDataType, tgtScans, baseFrame, smoothing);

% then run between scans motion compensation:
view = selectDataType(view, newDataType);
newnewDataType=['MotionComp_RefScan' num2str(baseScan)];
view = betweenScanMotComp(view, newnewDataType, baseScan, tgtScans);

disp('Finished all motion compensation. Final results in MotionComp_RefScan# data type.')

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [psnr, ssim] = compute_difference(im, imgt, SRF)
% COMPUTE_DIFFERENCE: compute image quality
%
% Input:
%     - im:   super-resolved image
%     - imgt: groundtruth high-resolution image
%     - SRF:  super-resolution factor
% Output:
%     - psnr: Peak signal-to-noise ratio
%     - ssim: Structural similarity index
%     - ifc:  Information fidelity criterion

% =========================================================================
% Retrieve only luminance channel
% =========================================================================
if size(im,3)>1
    im = rgb2ycbcr(im);
    im = im(:,:,1);
end

if size(imgt,3)>1
    imgt = rgb2ycbcr(imgt);
    imgt = imgt(:,:,1);
end

% =========================================================================
% Remove border pixels as some methods (e.g., A+) do not predict border pixels
% =========================================================================
cropPix     = SRF;
im          = shave(im, [cropPix, cropPix]);
imgt        = shave(imgt, [cropPix, cropPix]);

% Convert to double (with dynamic range 255)
im          = double(im);
imgt        = double(imgt);

% =========================================================================
% Compute Peak signal-to-noise ratio (PSNR)
% =========================================================================
mse = mean(mean((im - imgt).^2, 1), 2);
psnr = 10*log10(255*255/mse);

% =========================================================================
% Compute Structural similarity index (SSIM index)
% =========================================================================
[ssim, ~] = ssim_index(imgt, im);

% =========================================================================
% Compute information fidelity criterion (IFC)
% =========================================================================
% ifc = ifcvec(imgt, im);

end

function I = shave(I, border)
I = I(1+border(1):end-border(1), ...
      1+border(2):end-border(2), :, :);
end

author: tyshiwo


text:
function OutputFiles = bst_project_sources( ResultsFile, destSurfFile, isAbsoluteValues, isInteractive )
% BST_PROJECT_SOURCES: Project source files on a different surface (currents or timefreq).
%
% USAGE:  OutputFiles = bst_project_sources( ResultsFile, DestSurfFile, isAbsoluteValues=0, isInteractive=1 )
%
% INPUT:
%    - ResultsFile      : Relative path to sources file to reproject
%    - destSurfFile     : Relative path to destination surface file
%    - isAbsoluteValues : If 1, interpolate absolute values of the sources instead of relative values
%    - isInteractive    : If 1, displays questions and dialog boxes
%                         If 0, consider that it is running from the process interface

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

%% ===== PARSE INPUTS ======
if (nargin < 4) || isempty(isInteractive)
    isInteractive = 1;
end
if (nargin < 3) || isempty(isAbsoluteValues)
    isAbsoluteValues = 0;
end

%% ===== GROUP BY SURFACES =====
% Group the results files to process by the surface on which they were computed
% Objective: Computing only once the transformation for all the files in the same group
ResultsGroups = {};
SurfaceGroups = {};
nGroup = 0;
OutputFiles = {};
errMsg = [];
HeadModelType = [];
% Display progress bar
if isInteractive
    bst_progress('start', 'Project sources', 'Initialization...');
end
% Get file type: results or timefreq
isTimefreq = strcmpi(file_gettype(ResultsFile{1}), 'timefreq');
% For each sources file: get surface
for iRes = 1:length(ResultsFile)
    % Read surface file
    if isTimefreq
        ResMat = in_bst_timefreq(ResultsFile{iRes}, 0, 'SurfaceFile', 'DataType', 'DataFile', 'HeadModelType');
        % Check the data type: timefreq must be source/surface based, and no kernel-based file
        if ~strcmpi(ResMat.DataType, 'results')
            errMsg = 'Only cortical maps can be projected.';
            bst_report('Error', 'process_project_sources', ResultsFile{iRes}, errMsg);
            continue;
        elseif ~isempty(strfind(ResultsFile{iRes}, '_KERNEL_'))
            errMsg = 'Cannot re-project kernel-based time-frequency cortical maps.';
            bst_report('Error', 'process_project_sources', ResultsFile{iRes}, errMsg);
            continue;
        elseif isempty(ResMat.SurfaceFile) && ~isempty(ResMat.DataFile)
            ResAssocMat = in_bst_results(ResMat.DataFile, 0, 'SurfaceFile');
            ResMat.SurfaceFile = ResAssocMat.SurfaceFile;
        end
        % Load related source file
        if ~isempty(ResMat.DataFile) && isempty(ResMat.HeadModelType)
            [AssociateMat,AssociateFile] = in_bst_results(ResMat.DataFile, 0, 'HeadModelType');
            ResMat.HeadModelType = AssociateMat.HeadModelType;
        end
    else
        ResMat = in_bst_results(ResultsFile{iRes}, 0, 'SurfaceFile', 'HeadModelType');
    end

    % Default HeadModelType: surface
    if ~isfield(ResMat, 'HeadModelType') || isempty(ResMat.HeadModelType)
        ResMat.HeadModelType = 'surface';
    % Else : Check the type of grid (skip volume head models)
    elseif ismember(ResMat.HeadModelType, {'volume'})
        wrnMsg = ['To project source grids, see online tutorial "Group analysis: Subjects corergistration":' 10 ...
                  'https://neuroimage.usc.edu/brainstorm/Tutorials/CoregisterSubjects#Volume_source_models' 10 ...
                  'Skipping file: "' ResultsFile{iRes} '"...'];
        if isInteractive
            disp(wrnMsg);
        else
            bst_report('Error', 'process_project_sources', ResultsFile{iRes}, wrnMsg);
        end
        continue;
    end
    % Check head model type: must be the same for all the files
    if isempty(HeadModelType)
        HeadModelType = ResMat.HeadModelType;
    elseif ~isempty(ResMat.HeadModelType) && ~isequal(HeadModelType, ResMat.HeadModelType)
        errMsg = 'All the source files must be of the type (surface, volume or mixed).';
        bst_report('Error', 'process_project_sources', ResultsFile{iRes}, errMsg);
        continue;
    end

    % Associated surface not defined: error
    if isempty(ResMat.SurfaceFile)
        errMsg = 'Associated surface file is not defined.';
        bst_report('Error', 'process_project_sources', ResultsFile{iRes}, errMsg);
        continue;
    end
    % Check that it is not the destination surface
    if file_compare(ResMat.SurfaceFile, destSurfFile)
        if isInteractive
            disp(['BST> WARNING: Source and destination surfaces are the same for file: ' ResultsFile{iRes}]);
        else
            errMsg = 'Source and destination surfaces are the same.';
            bst_report('Error', 'process_project_sources', ResultsFile{iRes}, errMsg);
        end
        continue;
    end
    % Look for surface filename in SurfaceGroups
    iGroup = find(file_compare(ResMat.SurfaceFile, SurfaceGroups));
    % If does not exist yet: create it
    if isempty(iGroup)
        nGroup = nGroup + 1;
        SurfaceGroups{nGroup} = ResMat.SurfaceFile;
        ResultsGroups{nGroup} = ResultsFile(iRes);
    % Group exist: add file to the group
    else
        ResultsGroups{iGroup}{end+1} = ResultsFile{iRes};
    end
end
% If destination surface = source surface for all files
if (nGroup == 0)
    if isempty(errMsg)
        errMsg = ['Source and destination surfaces are the same for all the selected files.' 10 'Nothing to project...'];
    end
    if isInteractive
        bst_error(errMsg, 'Project sources', 0);
        bst_progress('stop');
    else
        bst_report('Error', 'process_project_sources', ResultsFile, errMsg);
    end
    return;
end
% Get protocol folders
ProtocolInfo = bst_get('ProtocolInfo');

%% ===== PROJECT SOURCES =====
isStopWarped = [];
iUpdatedStudies = [];
% Process each surface group
for iGroup = 1:nGroup
    % ===== GET INTERPOLATION =====
    srcSurfFile = SurfaceGroups{iGroup};
    if isInteractive
        bst_progress('start', 'Project sources', 'Loading surfaces...');
    end
    % Compute interpolation
    [WmatSurf, sSrcSubj, sDestSubj, srcSurfMat, destSurfMat, isStopWarped] = tess_interp_tess2tess(srcSurfFile, destSurfFile, isInteractive, isStopWarped);
    % Source subject and destination subject are the same
    isSameSubject = file_compare(sSrcSubj.FileName, sDestSubj.FileName);

    % ===== CREATE GROUP ANALYSIS SUBJECT =====
    % If src and dest subjects are not the same: create a "group analysis" subject
    if ~isSameSubject
        % If the destination is the group analysis subject: replace its name
        if strcmpi(sSrcSubj.Name, bst_get('DirDefaultSubject'))
            sSrcSubj.Name = bst_get('NormalizedSubjectName');
        end
        % If the destination is the group analysis subject: replace its name
        if strcmpi(sDestSubj.Name, bst_get('DirDefaultSubject'))
            sDestSubj = bst_get('NormalizedSubject');
            sDestSubj.Name = bst_get('NormalizedSubjectName');
        end
    end

    % ===== PROCESS EACH FILE =====
    nFile = length(ResultsGroups{iGroup});
    if isInteractive
        bst_progress('start', 'Project sources', 'Projecting sources...', 0, nFile);
    end
    % Process each results file in group
    for iFile = 1:nFile
        % Progress bar
        ResultsFile = ResultsGroups{iGroup}{iFile};
        if isInteractive
            bst_progress('inc', 1);
        end
        bst_progress('text', sprintf('Processing file #%d/%d: %s', iFile, nFile, ResultsFile));

        % ===== OUTPUT STUDY =====
        % Get source study
        [sSrcStudy, iSrcStudy] = bst_get('AnyFile', ResultsFile);
        % If result has to be save in "group analysis" subject
        if ~isSameSubject
            % New condition name
            NewCondition = strrep(sSrcStudy.Condition{1}, '@raw', '');
            % Get condition
            [sDestStudy, iDestStudy] = bst_get('StudyWithCondition', [sDestSubj.Name '/' NewCondition]);
            % Create condition if doesnt exist
            if isempty(iDestStudy)
                iDestStudy = db_add_condition(sDestSubj.Name, NewCondition, 0);
                if isempty(iDestStudy)
                    error(['Cannot create condition: "' normSubjName '/' NewCondition '".']);
                end
                sDestStudy = bst_get('Study', iDestStudy);
            end
        % Else: use the source study as output study
        else
            sDestStudy = sSrcStudy;
            iDestStudy = iSrcStudy;
        end

        % ===== LOAD INPUT FILE =====
        % Time-freq files
        if isTimefreq
            % Load file
            ResultsMat = in_bst_timefreq(ResultsFile, 0);
            ResFile = ResultsFile;
            % Change number of sources
            ResultsMat.RowNames = 1:size(ResultsMat.TF,1);
        % Source files: Read full
        else
            % Load file
            [ResultsMat,ResFile] = in_bst_results(ResultsFile, 1);
            % If it depends on the data file: try to use the comment of the data file
            if ~isempty(ResultsMat.DataFile) && strcmpi(file_gettype(ResultsMat.DataFile), 'data')
                % Find parent file in database
                [sStudy, iStudy, iData] = bst_get('DataFile', ResultsMat.DataFile);
                % If file was found: use its comment
                if ~isempty(sStudy) && ~isempty(sStudy.Data(iData).Comment)
                    ResultsMat.Comment = sStudy.Data(iData).Comment;
                end
            end
            % Unconstrained sources: Make a flat map
            if (ResultsMat.nComponents ~= 1)
                % Display warning before flattening, as this is not obvious
                strWarn = 'Unconstrained source maps are flattened (norm of the three orientations) before projection.';
                bst_report('Warning', 'process_project_sources', ResultsFile, strWarn);
                disp(['BST> Warning: ' strWarn]);
                % Apply norm of the three orientations
                ResultsMat = process_source_flat('Compute', ResultsMat);
            % Compute absolute values
            elseif isAbsoluteValues
                ResultsMat.ImageGridAmp = abs(ResultsMat.ImageGridAmp);
            end

            ResultsMat.ChannelFlag = [];
            ResultsMat.GoodChannel = [];
        end
        % Remove link with original file
        ResultsMat.DataFile = [];
        % Check if the file was reprojected on an atlas
        if isfield(ResultsMat, 'Atlas') && ~isempty(ResultsMat.Atlas)
            wrnMsg = ['Cannot process atlas-based source files: Skipping file "' ResultsFile '"...'];
            if isInteractive
                disp(wrnMsg);
            else
                bst_report('Error', 'process_project_sources', ResultsFile, wrnMsg);
            end
            continue;
        end

        % ===== INTERPOLATION MATRIX =====
        % Surface source model: Simply use the surface-surface interpolation
        if strcmpi(HeadModelType, 'surface')
            Wmat = WmatSurf;
        % Mixed source models: Compute volume grid interpolations
        elseif strcmpi(HeadModelType, 'mixed')
            % Load MRI files
            sMriSrc  = in_mri_bst(sSrcSubj.Anatomy(sSrcSubj.iAnatomy).FileName);
            sMriDest = in_mri_bst(sDestSubj.Anatomy(sDestSubj.iAnatomy).FileName);
            % Compute interpolation
            [Wmat, destGridAtlas, destGridLoc, destGridOrient] = tess_interp_mixed(ResultsMat, WmatSurf, srcSurfMat, destSurfMat, sMriSrc, sMriDest, isInteractive);
            % Update output structure
            ResultsMat.GridAtlas  = destGridAtlas;
            ResultsMat.GridLoc    = destGridLoc;
            ResultsMat.GridOrient = destGridOrient;
%             % Check if there is a "Source model" atlas available
%             iModelSrc  = find(strcmpi({srcSurfMat.Atlas.Name}, 'Source model'));
%             iModelDest = find(strcmpi({destSurfMat.Atlas.Name}, 'Source model'));
%             if isempty(iModelDest) && ~isempty(iModelSrc)
%                 destSurfMat.Atlas(end+1) = srcSurfMat.Atlas(iModelSrc);
%             end
        else
            error(['Unsupported head model type: ' HeadModelType]);
        end

        % ===== PROJECT SOURCE MAPS =====
        % Time-freq file
        if isTimefreq
            % Apply interpolation matrix
            tmpTF = zeros(size(Wmat,1), size(ResultsMat.TF,2), size(ResultsMat.TF,3));
            for iFreq = 1:size(ResultsMat.TF,3)
                tmpTF(:,:,iFreq) = Wmat * ResultsMat.TF(:,:,iFreq);
            end
            ResultsMat.TF = tmpTF;
            % PAC: Apply interpolation to all measures
            if isfield(ResultsMat, 'sPAC') && ~isempty(ResultsMat.sPAC)
                if isfield(ResultsMat.sPAC, 'NestingFreq') && ~isempty(ResultsMat.sPAC.NestingFreq)
                    ResultsMat.sPAC.NestingFreq = Wmat * ResultsMat.sPAC.NestingFreq;
                end
                if isfield(ResultsMat.sPAC, 'NestedFreq') && ~isempty(ResultsMat.sPAC.NestedFreq)
                    ResultsMat.sPAC.NestedFreq = Wmat * ResultsMat.sPAC.NestedFreq;
                end
                if isfield(ResultsMat.sPAC, 'DirectPAC') && ~isempty(ResultsMat.sPAC.DirectPAC)
                    tmpTF = zeros(size(Wmat,1), size(ResultsMat.sPAC.DirectPAC,2), size(ResultsMat.sPAC.DirectPAC,3), size(ResultsMat.sPAC.DirectPAC,4));
                    for iLow = 1:size(ResultsMat.sPAC.DirectPAC,3)
                        for iHigh = 1:size(ResultsMat.sPAC.DirectPAC,4)
                            tmpTF(:,:,iLow,iHigh) = Wmat * ResultsMat.sPAC.DirectPAC(:,:,iLow,iHigh);
                        end
                    end
                    ResultsMat.sPAC.DirectPAC = tmpTF;
                end
            end
            ResultsMat.HeadModelType = HeadModelType;
            if ~isempty(ResultsMat.GridLoc) && (length(ResultsMat.RowNames) ~= size(ResultsMat.GridLoc,1))
                ResultsMat.RowNames = 1:size(ResultsMat.GridLoc,1);
            elseif isequal(ResultsMat.DataType, 'results') && isnumeric(ResultsMat.RowNames) && (length(ResultsMat.RowNames) ~= size(ResultsMat.TF,1))
                ResultsMat.RowNames = 1:size(ResultsMat.TF,1);
            end
        % Source file
        else
            % Apply interpolation matrix
            ResultsMat.ImageGridAmp = muliplyInterp(Wmat, double(ResultsMat.ImageGridAmp), ResultsMat.nComponents);

            % Apply interpolation to standart deviation matrix
            if isfield(ResultsMat, 'Std') && ~isempty(ResultsMat.Std)
                ResultsMat.Std = muliplyInterp(Wmat, double(ResultsMat.Std), ResultsMat.nComponents);
            end
            ResultsMat.ImagingKernel = [];
        end

        % === SAVE NEW RESULTS ===
        % Get source filename
        [tmp__, oldBaseName] = bst_fileparts(ResFile);
        % Remove KERNEL tag for saving full source files
        if ~isempty(strfind(oldBaseName, '_KERNEL_')) && isfield(ResultsMat, 'ImageGridAmp') && ~isempty(ResultsMat.ImageGridAmp)
            oldBaseName = strrep(oldBaseName, '_KERNEL', '');
        end
        % Prepare structure to be saved
        ResultsMat.SurfaceFile = destSurfFile;
        if ~isSameSubject
            ResultsMat.Comment = [sSrcSubj.Name '/' ResultsMat.Comment];
            newResultsFile = sprintf('%s_%s.mat', oldBaseName, file_standardize(sSrcSubj.Name));
        else
            ResultsMat.Comment = [ResultsMat.Comment ' | ' destSurfMat.Comment];
            newResultsFile = sprintf('%s_%dV.mat', oldBaseName, length(destSurfMat.Vertices));
        end
        % Surface file
        ResultsMat.SurfaceFile = file_win2unix(destSurfFile);
        % History: project source
        ResultsMat = bst_history('add', ResultsMat, 'project', ['Project sources: ' srcSurfFile ' => ' ResultsMat.SurfaceFile]);
        % Build full filename
        newResultsFileFull = bst_fullfile(ProtocolInfo.STUDIES, bst_fileparts(sDestStudy.FileName), newResultsFile);
        newResultsFileFull = file_unique(newResultsFileFull);
        newResultsFile = file_short(newResultsFileFull);
        % Save new results file
        bst_save(newResultsFileFull, ResultsMat, 'v6');

        % === ADD FILE IN DATABASE ===
        % Create Results/Timefreq structure for database
        if isTimefreq
            sNewResults = db_template('Timefreq');
            sNewResults.FileName = newResultsFile;
            sNewResults.Comment  = ResultsMat.Comment;
            sNewResults.DataFile = ResultsMat.DataFile;
            sNewResults.DataType = ResultsMat.DataType;
            % If filename already exists in this study
            iExistingRes = find(file_compare({sDestStudy.Timefreq.FileName}, newResultsFile));
            if ~isempty(iExistingRes)
                % Replace previous Results
                sDestStudy.Timefreq(iExistingRes) = sNewResults;
            else
                % Add new result
                sDestStudy.Timefreq(end + 1) = sNewResults;
            end
        else
            sNewResults = db_template('Results');
            sNewResults.FileName = newResultsFile;
            sNewResults.Comment  = ResultsMat.Comment;
            sNewResults.DataFile = ResultsMat.DataFile;
            sNewResults.isLink   = 0;
            sNewResults.HeadModelType = HeadModelType;
            % If filename already exists in this study
            iExistingRes = find(file_compare({sDestStudy.Result.FileName}, newResultsFile));
            if ~isempty(iExistingRes)
                % Replace previous Results
                sDestStudy.Result(iExistingRes) = sNewResults;
            else
                % Add new result
                sDestStudy.Result(end + 1) = sNewResults;
            end
        end
        % Update study in database
        bst_set('Study', iDestStudy, sDestStudy);
        iUpdatedStudies = [iUpdatedStudies, iDestStudy];
        % Add to list of returned files
        OutputFiles{end+1} = newResultsFile;
    end
end

%% ===== UDPATE DISPLAY =====
if isInteractive
    bst_progress('stop');
end
if isempty(OutputFiles)
    return;
end
% Update node
panel_protocols('UpdateNode', 'Study', unique(iUpdatedStudies));
% Select first output study
panel_protocols('SelectStudyNode', iUpdatedStudies(1));
% Select first output file
panel_protocols('SelectNode', [], OutputFiles{1});
% Save database
db_save();

end

%% ===== APPLY INTERPOLATION MATRIX =====
function B = muliplyInterp(W, A, nComponents)
    switch (nComponents)
        case {0, 1}
            B = double(W * A);
        case 2
            B = zeros(2 * size(W,1), size(A,2));
            B(1:2:end,:) = double(W * A(1:2:end,:));
            B(2:2:end,:) = double(W * A(2:2:end,:));
        case 3
            B = zeros(3 * size(W,1), size(A,2));
            B(1:3:end,:) = double(W * A(1:3:end,:));
            B(2:3:end,:) = double(W * A(2:3:end,:));
            B(3:3:end,:) = double(W * A(3:3:end,:));
    end
    B = double(B);
end

author: brainstorm-tools


text:
function sv_startSchuster(mCatalog, hFigure, bMap, rContainer)
% function sv_startSchuster(mCatalog, hFigure, bMap, rContainer)
% ------------------------------------------------------
% Starts seismicity variation analysis toolbox: Schuster Test on Mc
%
% Input parameters:
%   mCatalog      Earthquake catalog to be analyzed
%   hFigure       Handle of figure where the user should select the grid (i.e. the seismicity map)
%   bMap          Map/cross-section switch. If the testing is carried out on a map set bMap = 1,
%                 on a cross-section set bMap = 0
%   rContainer    Just a container (structure proposed) to store any variables
%
% J. Woessner; woessner@seismo.ifg.ethz.ch
% April 4, 2003

global bDebug;
if bDebug
    report_this_filefun(mfilename('fullpath'));
end

% Launch GUI
hMenuFig = sv_gridSchuster(1);

% Set up parameter struct and store input parameters
params.mCatalog = mCatalog;
params.bMap = bMap;
params.rContainer = rContainer;

% Analyze Output
if ~ishandle(hMenuFig)
    answer = 0;
else
    handles = guidata(hMenuFig);
    answer = handles.answer;
    % OK pressed
    if answer == 1
        % Get the values from figure
        %params.bNumber = get(handles.radNumber, 'Value');   % deprecated
        if get(handles.radNumber, 'Value') == 1
            params.nGriddingMode = 0;   % Constant number
        elseif get(handles.radRadius, 'Value') == 1
            params.nGriddingMode = 1;   % Constant radius
        else
            params.nGriddingMode = 2;   % Rectangle mode
        end
        % Setting up the params struct array with variables
        params.nNumberEvents = str2double(get(handles.txtNumber, 'String'));
        params.fMaxRadius = str2double(get(handles.txtMaxRadius, 'String'));
        params.fRadius = str2double(get(handles.txtRadius, 'String'));
        params.bGridEntireArea = get(handles.chkGridEntireArea, 'Value');
        params.fSpacingHorizontal = str2double(get(handles.txtSpacingHorizontal, 'String')); % Grid spacing variable
        params.fSpacingDepth = str2double(get(handles.txtSpacingDepth, 'String'));           % Grid spacing variable
        params.fSizeRectHorizontal = str2double(get(handles.txtSizeRectHorizontal, 'String')); % Rectangular selection instead radius params.fRadius
        params.fSizeRectDepth = str2double(get(handles.txtSizeRectDepth, 'String'));           % Rectangular selection instead radius params.fRadius
        params.nMinimumNumber = str2double(get(handles.txtMinimumNumber, 'String'));
        params.bTimePeriod = get(handles.chkTimePeriod,'Value');
        params.fTimePeriod = str2double(get(handles.txtTimePeriod,'String')); % Period lengths to be compared in days
        params.bTstart = get(handles.chkTstart,'Value'); % Check for starting time of temporal mapping
        params.fTstart = str2double(get(handles.txtTstart,'String')); % Starting time for temporal mapping
        params.bBstnum = get(handles.chkBstnum,'Value'); % Check for boostrap sampling
        params.fBstnum = str2double(get(handles.txtBstnum,'String')); % Number of bootstrap samples
        params.fBinning = str2double(get(handles.txtBinsize,'String')); % Bin size for magnitude binning
        params.fStartMag = str2double(get(handles.txtStartMag,'String')); % Lower magnitude to perform Schuster test
        params.fEndMag = str2double(get(handles.txtEndMag,'String')); % Last magnitude to perform Schuster test
        params.sComment = get(handles.txtComment, 'String'); % Comment on calculation

        %%%%%%%%%%%%%%%%%%%%%%%%%% Check variable settings %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %
        %  if ~params.bMinMagMc
        %      params.fMinMag = str2double(get(handles.txtMinMag, 'String'));
        %    else
        %      params.fMinMag = 1;
        %    end
        %    params.fMaxMag = str2double(get(handles.txtMaxMag, 'String'));
        bSaveParameter = get(handles.chkSaveParameter, 'Value');
        if bSaveParameter
            sSaveParameter = get(handles.lblSaveParameter, 'String');
        end

        % Remove figure from memory
        delete(hMenuFig);

        % Select grid
        [params.mPolygon, params.vX, params.vY, params.vUsedNodes] = ex_selectgrid(hFigure, params.fSpacingHorizontal, params.fSpacingDepth, params.bGridEntireArea);

        % Validate polygonsize
        if length(params.vX) < 4  ||  length(params.vY) < 4
            errordlg('Selection is too small. Please select a larger polygon.');
            return;
        end

        % Add parameter to params.sComment
        params.sComment = [params.sComment ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.'...
                    ', Time period ' num2str(params.fTimePeriod) ' d'];
        if bSaveParameter
            % Save paramters
            save(sSaveParameter, 'params');
        else
            % General calculations

            % Create Indices to catalog
            [params.caNodeIndices] = ex_CreateIndexCatalog(params.mCatalog, params.mPolygon, params.bMap, params.nGriddingMode, ...
                params.nNumberEvents, params.fRadius, params.fSizeRectHorizontal, params.fSizeRectDepth);

            %%%%% Temporarily calculations for %%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Perform the calculation
            [params] = sv_calcSchuster(params);

            sv_result(params);
        end
    else
        delete(hMenuFig);
    end
end


text:
classdef PTKSettings < CoreBaseClass
    % PTKSettings. Part of the internal gui for the TD MIM Toolkit.
    %
    %     You should not use this class within your own code. It is intended to
    %     be used internally within the gui of the TD MIM Toolkit.
    %
    %     PTKSettings stores gui application settings between sessions. This
    %     class is stored on disk so that settings persist betweek sessions.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        OverlayOpacity = 50
        ImageInfo = []
        SliceNumber = [1 1 1]
        SaveImagePath = ''
        ShowTextLabels = true
        CurrentMarkerColour = 1
        ScreenPosition
        PatientBrowserScreenPosition
        DeveloperMode = false
        LastUidForPatientMap
        LastMarkerSetForPatientMap
    end

    properties (Transient, Access = private)
        SettingsFilename
    end

    methods (Static)
        function settings = LoadSettings(app_def, reporting)
            try
                settings_filename = app_def.GetSettingsFilePath;
                if exist(settings_filename, 'file')
                    settings_struct = MimDiskUtilities.Load(settings_filename);
                    settings = settings_struct.settings;
                else
                    reporting.Log('No settings file found. Will create new one on exit');
                    settings = PTKSettings;
                end
                settings.SettingsFilename = settings_filename;
                if isempty(settings.LastMarkerSetForPatientMap)
                    settings.LastMarkerSetForPatientMap = containers.Map;
                end
                if isempty(settings.LastUidForPatientMap)
                    settings.LastUidForPatientMap = containers.Map;
                end

            catch ex
                reporting.ErrorFromException('PTKSettings:FailedtoLoadSettingsFile', ['Error when loading settings file ' settings_filename '. Try deleting this file.'], ex);
            end

        end

    end

    methods
        function obj = PTKSettings
            obj.LastUidForPatientMap = containers.Map;
            obj.LastMarkerSetForPatientMap = containers.Map;
        end

        function ApplySettingsToGui(obj, gui, viewer_panel)
            gui.DeveloperMode = obj.DeveloperMode;
            viewer_panel.OverlayOpacity = obj.OverlayOpacity;
            viewer_panel.SliceNumber = obj.SliceNumber;
            viewer_panel.MarkerImageDisplayParameters.ShowLabels = obj.ShowTextLabels;
            viewer_panel.NewMarkerColour = obj.CurrentMarkerColour;
        end

        function UpdateSettingsFromGui(obj, gui, viewer_panel)
            obj.DeveloperMode = gui.DeveloperMode;
            obj.OverlayOpacity = viewer_panel.OverlayOpacity;
            obj.SliceNumber = viewer_panel.SliceNumber;
            obj.ShowTextLabels = viewer_panel.MarkerImageDisplayParameters.ShowLabels;
            obj.CurrentMarkerColour = viewer_panel.NewMarkerColour;
        end

        function SetLastImageInfo(obj, image_info, reporting)
            if ~isequal(image_info, obj.ImageInfo)
                obj.ImageInfo = image_info;
                obj.SaveSettings(reporting);
            end
        end

        function AddLastPatientUid(obj, patient_id, series_uid)
            if isempty(obj.LastUidForPatientMap)
                obj.LastUidForPatientMap = containers.Map;
            end
            obj.LastUidForPatientMap(patient_id) = series_uid;
        end

        function AddLastMarkerSet(obj, series_uid, marker_set_name)
            if isempty(obj.LastMarkerSetForPatientMap)
                obj.LastMarkerSetForPatientMap = containers.Map;
            end
            obj.LastMarkerSetForPatientMap(series_uid) = marker_set_name;
        end

        function RemoveLastMarkerSet(obj, series_uid, marker_set_name)
            if ~isempty(obj.LastMarkerSetForPatientMap) && obj.LastMarkerSetForPatientMap.isKey(series_uid) && strcmp(obj.LastMarkerSetForPatientMap(series_uid), marker_set_name)
                obj.LastMarkerSetForPatientMap.remove(series_uid);
            end
        end

        function series_uid = GetLastPatientUid(obj, patient_id)
            if obj.LastUidForPatientMap.isKey(patient_id)
                series_uid = obj.LastUidForPatientMap(patient_id);
            else
                series_uid = [];
            end
        end

        function marker_set_name = GetLastMarkerSetName(obj, series_uid)
            if isempty(obj.LastMarkerSetForPatientMap)
                obj.LastMarkerSetForPatientMap = containers.Map;
            end
            if obj.LastMarkerSetForPatientMap.isKey(series_uid)
                marker_set_name = obj.LastMarkerSetForPatientMap(series_uid);
            else
                marker_set_name = [];
            end
        end

        function RemoveLastPatientUid(obj, series_uids)
            for key = obj.LastUidForPatientMap.keys
                if any(strcmp(obj.LastUidForPatientMap(key{1}), series_uids))
                    obj.LastUidForPatientMap.remove(key{1});
                end
            end
            for key = obj.LastMarkerSetForPatientMap.keys
                if any(strcmp(obj.LastMarkerSetForPatientMap(key{1}), series_uids))
                    obj.LastMarkerSetForPatientMap.remove(key{1});
                end
            end

            if ~isempty(obj.ImageInfo) && any(strcmp(obj.ImageInfo.ImageUid, series_uids))
                obj.ImageInfo = [];
            end
        end

        function SetLastSaveImagePath(obj, image_path, reporting)
            if ~isempty(image_path) && ischar(image_path) && ~strcmp(image_path, obj.SaveImagePath)
                obj.SaveImagePath = image_path;
                obj.SaveSettings(reporting);
            end
        end

        function SetPosition(obj, main_screen_position, patient_browser_position)
            obj.ScreenPosition = main_screen_position;
            obj.PatientBrowserScreenPosition = patient_browser_position;
        end

        function SaveSettings(obj, reporting)

            settings_filename = obj.SettingsFilename;

            try
                value = [];
                value.settings = obj;
                if isempty(settings_filename)
                    reporting.ShowWarning('PTKSettings:NoSettingsFilename', 'The settings file could not be saved as the settings filename was not known.', []);
                else
                    MimDiskUtilities.Save(settings_filename, value);
                end
            catch ex
                reporting.ErrorFromException('PTKSettings:FailedtoSaveSettingsFile', ['Unable to save settings file ' settings_filename], ex);
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [mph] = mach2mph(mach)
% Convert speed from mach (at standard temp and pressure!) to miles per hour.
% Chad A. Greene 2012
mph = mach*767.2691481747;

author: Sable

repo: mcbench-benchmarks


text:
function [A,h] = build_A_h_3d(dphi,cert,phiGrad,transformatioinModel)
% BUILD_A_H_LINEAR3D Builds the equation system A*p = h
%
% INPUT ARGUMENTS
% dphi                  - Phase-difference
% cert                  - Certainty
% phiGrad               - Phase gradient
% transformationModel   - Transformation model (translation or affine)
%
% OPTIONAL INPUT ARGUMENTS
% N/A
%
% OUTPUT ARGUMENTS
% A 					- A matrix
% h 					- h vector
%
% See "Phase-Based Multidimensional Volume Registration" by Hemmendorf et al
% or "PHASE BASED VOLUME REGISTRATION USING CUDA" by Eklund et al for a
% detailed description of how the equation system is set up. Note that here
% we use a differenct S(x)p then the one used in the papers. This is done
% in order to be consistent with functions build_G_h_linear2d and
% build_G_h_linear3d.

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

sz = size(dphi(:,:,:,1));

dphiX = vec(dphi(:,:,:,1));
dphiY = vec(dphi(:,:,:,2));
dphiZ = vec(dphi(:,:,:,3));
certX = vec(cert(:,:,:,1));
certY = vec(cert(:,:,:,2));
certZ = vec(cert(:,:,:,3));
phiGradX = vec(phiGrad(:,:,:,1));
phiGradY = vec(phiGrad(:,:,:,2));
phiGradZ = vec(phiGrad(:,:,:,3));

switch transformatioinModel
    case 'translation'
        A = [   sum(certX.*phiGradX.^2),    0,                          0; ...
                0,                          sum(certY.*phiGradY.^2),    0; ...
                0,                          0,                          sum(certZ.*phiGradZ.^2)];
        h = [   sum(certX.*dphiX.*phiGradX); ...
                sum(certY.*dphiY.*phiGradY); ...
                sum(certZ.*dphiZ.*phiGradZ)];
    case {'rigid','affine'}
        [x,y,z] = meshgrid(1:sz(2),1:sz(1),1:sz(3));
        x = x - sz(2)/2 - 0.5;
        y = y - sz(1)/2 - 0.5;
        z = z - sz(3)/2 - 0.5;
        x = x(:);
        y = y(:);
        z = z(:);

        A = [   sum(certX.*phiGradX.^2.*x.^2),  sum(certX.*phiGradX.^2.*x.*y),  sum(certX.*phiGradX.^2.*x.*z)       0,                                  0,                              0,                              0,                              0,                              0,                              sum(certX.*phiGradX.^2.*x), 0,                          0; ...
                0,                              sum(certX.*phiGradX.^2.*y.^2),  sum(certX.*phiGradX.^2.*y.*z)       0,                                  0,                              0,                              0,                              0,                              0,                              sum(certX.*phiGradX.^2.*y), 0,                          0; ...
                0,                              0,                              sum(certX.*phiGradX.^2.*z.^2)       0,                                  0,                              0,                              0,                              0,                              0,                              sum(certX.*phiGradX.^2.*z), 0,                          0; ...
                0,                              0,                              0,                                  sum(certY.*phiGradY.^2.*x.^2),      sum(certY.*phiGradY.^2.*x.*y),  sum(certY.*phiGradY.^2.*x.*z),  0,                              0,                              0,                              0,                          sum(certY.*phiGradY.^2.*x), 0; ...
                0,                              0,                              0,                                  0,                                  sum(certY.*phiGradY.^2.*y.^2),  sum(certY.*phiGradY.^2.*y.*z),  0,                              0,                              0,                              0,                          sum(certY.*phiGradY.^2.*y), 0; ...
                0,                              0,                              0,                                  0,                                  0,                              sum(certY.*phiGradY.^2.*z.^2),  0,                              0,                              0,                              0,                          sum(certY.*phiGradY.^2.*z), 0; ...
                0,                              0,                              0,                                  0,                                  0,                              0,                              sum(certZ.*phiGradZ.^2.*x.^2),  sum(certZ.*phiGradZ.^2.*x.*y),  sum(certZ.*phiGradZ.^2.*x.*z),  0,                          0,                          sum(certZ.*phiGradZ.^2.*x); ...
                0,                              0,                              0,                                  0,                                  0,                              0,                              0,                              sum(certZ.*phiGradZ.^2.*y.^2),  sum(certZ.*phiGradZ.^2.*y.*z),  0,                          0,                          sum(certZ.*phiGradZ.^2.*y); ...
                0,                              0,                              0,                                  0,                                  0,                              0,                              0,                              0,                              sum(certZ.*phiGradZ.^2.*z.^2),  0,                          0,                          sum(certZ.*phiGradZ.^2.*z); ...
                0,                              0,                              0,                                  0,                                  0,                              0,                              0,                              0,                              0,                              sum(certX.*phiGradX.^2),    0,                          0; ...
                0,                              0,                              0,                                  0,                                  0,                              0,                              0,                              0,                              0,                              0,                          sum(certY.*phiGradY.^2),    0; ...
                0,                              0,                              0,                                  0,                                  0,                              0,                              0,                              0,                              0,                              0,                          0,                          sum(certZ.*phiGradZ.^2)];

        h = [   sum(certX.*dphiX.*phiGradX.*x); ...
                sum(certX.*dphiX.*phiGradX.*y); ...
                sum(certX.*dphiX.*phiGradX.*z); ...
                sum(certY.*dphiY.*phiGradY.*x); ...
                sum(certY.*dphiY.*phiGradY.*y); ...
                sum(certY.*dphiY.*phiGradY.*z); ...
                sum(certZ.*dphiZ.*phiGradZ.*x); ...
                sum(certZ.*dphiZ.*phiGradZ.*y); ...
                sum(certZ.*dphiZ.*phiGradZ.*z); ...
                sum(certX.*dphiX.*phiGradX); ...
                sum(certY.*dphiY.*phiGradY); ...
                sum(certZ.*dphiZ.*phiGradZ)];

        A = A + transpose(A) - diag(diag(A));
end

author: fordanic

repo: image-registration


text:
function varargout = exp(varargin)
%EXP (overloaded)

switch class(varargin{1})

    case 'double' % What is the numerical value of this argument (needed for displays etc)
        % SHOULD NEVER HAPPEN, THIS SHOULD BE CAUGHT BY BUILT-IN
        error('Overloaded SDPVAR/NORM CALLED WITH DOUBLE. Report error')

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.
        if length(varargin{1}) == 1
            varargout{1} = yalmip('addEvalVariable',mfilename,varargin{1});
        else
            y = [];
            for i = 1:length(varargin{1})
                y = [y;yalmip('addEvalVariable',mfilename,extsubsref(varargin{1},i))];
            end
            varargout{1} = y;
        end

    case 'char' % YALMIP sends 'model' when it wants the epigraph or hypograph
        switch varargin{1}
            case 'graph'
                t = varargin{2};
                X = varargin{3};

                % This is different from so called extended operators
                % Just do it!
                F = SetupEvaluationVariable(varargin{:});

                % Now add your own code, such as domain constraints.
                % Exponential does not need any domain constraint.

                % Let YALMIP know about convexity etc
                varargout{1} = F;
                varargout{2} = struct('convexity','convex','monotonicity','increasing','definiteness','positive');
                varargout{3} = X;

            case 'milp'
                    varargout{1} = [];
                    varargout{2} = [];
                    varargout{3} = [];
            otherwise
                error('SDPVAR/EXP called with CHAR argument?');
        end
    otherwise
        error('SDPVAR/EXP called with CHAR argument?');
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [uSol, NuSol] = startMultistep(K, dt, L, Nc, Nv, ~, S, uInit, NuInit)
%STARTMULTISTEP   Get enough initial data when using a multistep scheme.
%    [USOL, NUSOL] = STARTMULTISTEP(K, dt, L, NC, NV, pref, S, uInit, NuInit)
%    uses a one-step algorithm to start a multistep scheme.

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

% Note for developers: multistep algorithms are started with the one-step LIRK4
% time-stepping scheme.

% Set-up:
q = K.steps;                      % number of steps
nVars = S.numVars;                % number of unknown functions
%N = sqrt(size(uSol{1}, 1)/nVars); % number of grid points

% Create a cell-array to store the coefficients at the Q steps:
uSol = cell(q, 1);
NuSol = cell(q, 1);

% Store the initial conidition in the last column:
uSol{q} = uInit{1};
NuSol{q} = NuInit{1};

% Set-up the scheme:
K = imex('lirk4');
schemeCoeffs = computeCoeffs(K, dt, L, [], S);

% Do (Q-1) steps of LIRK4:
uOld = uInit;
NuOld = NuInit;
for j = 1:(q-1)
    [uNew, NuNew] = oneStep(K, dt, schemeCoeffs, Nc, Nv, nVars, S, uOld, NuOld);
    uSol{q-j} = uNew{1};
    NuSol{q-j} = NuNew{1};
    uOld = uNew;
    NuOld = NuNew;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%%% movobj = convert_3dresults2mov(3dmatrix,3dmatrix,3dmatrix,3dmatrix,struct)
% I - Input sequence
% L - Low-rank sequence
% S - Sparse sequence
% O - Outliers sequence
%
function [movobj] = convert_3dresults2mov(I,L,S,O,nframes)
  %warning('off','all');
  displog('Converting results to movie...');

  if(~isempty(I)) movobj.I(1:nframes) = struct('cdata', [], 'colormap', []); else movobj.I = []; end
  if(~isempty(L)) movobj.L(1:nframes) = struct('cdata', [], 'colormap', []); else movobj.L = []; end
  if(~isempty(S)) movobj.S(1:nframes) = struct('cdata', [], 'colormap', []); else movobj.S = []; end
  if(~isempty(O)) movobj.O(1:nframes) = struct('cdata', [], 'colormap', []); else movobj.O = []; end

  for i = 1 : nframes
    %
    % Convert input
    %
    if(~isempty(I))
      Input = I(:,:,i);
      Input = mat2gray(Input);
      Input = im2uint8(Input);
      rgbImage = repmat(Input,[1 1 3]);
      movobj.I(i).cdata = rgbImage;
    end
    %
    % Convert low-rank results
    %
    if(~isempty(L))
      LowRank = L(:,:,i);
      LowRank = mat2gray(LowRank);
      LowRank = im2uint8(LowRank);
      rgbImage = repmat(LowRank,[1 1 3]);
      movobj.L(i).cdata = rgbImage;
    end
    %
    % Store sparse results
    %
    if(~isempty(S))
      Sparse = S(:,:,i);
      Sparse = mat2gray(Sparse);
      Sparse = im2uint8(Sparse);
      rgbImage = repmat(Sparse,[1 1 3]);
      movobj.S(i).cdata = rgbImage;
    end
    %
    % Store foreground results
    %
    if(~isempty(O))
      Outlier = O(:,:,i);
      Outlier = mat2gray(Outlier);
      Outlier = im2uint8(Outlier);
      Outlier = medfilt2(Outlier, [5 5]);
      rgbImage = repmat(Outlier,[1 1 3]);
      movobj.O(i).cdata = rgbImage;
    end
  end

  %warning('on','all');
end

author: andrewssobral

repo: lrslibrary


text:
% timtopo()   - plot all channels of a data epoch on the same axis
%               and map its scalp map(s) at selected latencies.
% Usage:
%  >> timtopo(data, chan_locs);
%  >> timtopo(data, chan_locs, 'key', 'val', ...);
% Inputs:
%  data       = (channels,frames) single-epoch data matrix
%  chan_locs  = channel location file or EEG.chanlocs structure.
%               See >> topoplot example for file format.
%
% Optional ordered inputs:
%  'limits'    = [minms maxms minval maxval] data limits for latency (in ms) and y-values
%                 (assumes uV) {default|0 -> use [0 npts-1 data_min data_max];
%                 else [minms maxms] or [minms maxms 0 0] -> use
%                [minms maxms data_min data_max]
%  'plottimes' = [vector] latencies (in ms) at which to plot scalp maps
%                {default|NaN -> latency of maximum variance}
% 'title'      = [string] plot title {default|0 -> none}
% 'plotchans'  = vector of data channel(s) to plot. Note that this does not
%                affect scalp topographies {default|0 -> all}
% 'voffsets'   = vector of (plotting-unit) distances vertical lines should extend
%                above the data (in special cases) {default -> all = standard}
%
% Optional keyword, arg pair inputs (must come after the above):
% 'topokey','val' = optional topoplot() scalp map plotting arguments. See >> help topoplot
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 1-10-98
%
% See also: envtopo(), topoplot()

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

% 5-31-00 added o-time line and possibility of plotting 1 channel -sm & mw
% 11-02-99 added maplimits arg -sm
% 01-22-01 added to help message -sm
% 01-25-02 reformated help & license, added link -ad
% 03-15-02 add all topoplot options -ad

function M = timtopo(data, chan_locs, varargin)

MAX_TOPOS = 24;

if nargin < 1 %should this be 2?
   help timtopo;
   return
end

[chans,frames] = size(data);
icadefs;

%
% if nargin > 2
%     if ischar(limits)
%         varargin = { limits,plottimes,titl,plotchans,voffsets varargin{:} };
%
%     else
%         varargin = { 'limits' limits 'plottimes' plottimes 'title' titl 'plotchans' plotchans 'voffsets' voffsets varargin{:} };
%     end;
% end;

if nargin > 2 && ~ischar(varargin{1})
   options = {};
   if length(varargin) > 0, options = { options{:} 'limits' varargin{1} }; end;
   if length(varargin) > 1, options = { options{:} 'plottimes' varargin{2} }; end;
   if length(varargin) > 2, options = { options{:} 'title'      varargin{3} }; end;
   if length(varargin) > 3, options = { options{:} 'plotchans' varargin{4} }; end;
   if length(varargin) > 4, options = { options{:} 'voffsets'     varargin{5} }; end;
   if length(varargin) > 5, options = { options{:} varargin{6:end} }; end;
else
   options = varargin;
end;

fieldlist = { 'limits'        'real'     []                       0;
              'plottimes'     'real'     []                       [];
              'title'         'string'   []                       '';
              'plotchans'     'integer'  [1:size(data,1)]         0;
              'voffsets'      'real'     []                       [] ;};
[g topoargs] = finputcheck(options, fieldlist, 'timtopo', 'ignore');

if ischar(g), error(g); end;
%Set Defaults
if isempty(g.title), g.title = ''; end;
if isempty(g.voffsets) || g.voffsets == 0, g.voffsets = zeros(1,MAX_TOPOS); end;
if isempty(g.plotchans) || isequal(g.plotchans,0), g.plotchans = 1:chans; end;
plottimes_set=1;   % flag variable
if isempty(g.plottimes) || any(isnan(g.plottimes)), plottimes_set = 0;end;
limitset = 0; %flag variable
if isempty(g.limits), g.limits = 0; end;
if length(g.limits)>1, limitset = 1; end;

% if nargin < 7 | voffsets == 0
%   voffsets = zeros(1,MAX_TOPOS);
% end
%
% if nargin < 6
%    plotchans = 0;
% end
%
% if plotchans==0
%    plotchans = 1:chans;
% end
%
% if nargin < 5,
%    titl = '';     % DEFAULT NO TITLE
% end
%
% plottimes_set=1;   % flag variable
% if nargin< 4 | isempty(plottimes) | any(isnan(plottimes))
%    plottimes_set = 0;
% end
%
% limitset = 0;
% if nargin < 3,
%     limits = 0;
% elseif length(limits)>1
%     limitset = 1;
% end

if nargin < 2 %if first if-statement is changed to 2 should this be 3?
    chan_locs = 'chan.locs';  % DEFAULT CHAN_FILE
end
if isnumeric(chan_locs) && chan_locs == 0,
    chan_locs = 'chan.locs';  % DEFAULT CHAN_FILE
end

  %
  %%%%%%%%%%%%%%%%%%%%%%% Read and adjust limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % defaults: limits == 0 or [0 0 0 0]
  if ( length(g.limits) == 1 & g.limits==0) | (length(g.limits)==4 & ~any(g.limits))
    xmin=0;
    xmax=frames-1;
    ymin=min(min(data));
    ymax=max(max(data));
  elseif length(g.limits) == 2  % [minms maxms] only
    ymin=min(min(data));
    ymax=max(max(data));
    xmin = g.limits(1);
    xmax = g.limits(2);
 elseif length(g.limits) == 4
    xmin = g.limits(1);
    xmax = g.limits(2);
    if any(g.limits([3 4]))
      ymin = g.limits(3);
      ymax = g.limits(4);
    else % both 0
      ymin=min(min(data));
      ymax=max(max(data));
    end
  else
    fprintf('timtopo(): limits format not correct. See >> help timtopo.\n');
    return
  end;

  if xmax == 0 & xmin == 0,
    x = (0:1:frames-1);
    xmin = 0;
    xmax = frames-1;
  else
    dx = (xmax-xmin)/(frames-1);
    x=xmin*ones(1,frames)+dx*(0:frames-1); % compute x-values
  end;
  if xmax<=xmin,
      fprintf('timtopo() - in limits, maxms must be > minms.\n')
      return
  end

  if ymax == 0 & ymin == 0,
      ymax=max(max(data));
      ymin=min(min(data));
  end
  if ymax<=ymin,
      fprintf('timtopo() - in limits, maxval must be > minmval.\n')
      return
  end

sampint = (xmax-xmin)/(frames-1); % sampling interval = 1000/srate;
x = xmin:sampint:xmax;   % make vector of x-values

%
%%%%%%%%%%%%%%% Compute plot times/frames %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if plottimes_set == 0
  [mx plotframes] = max(sum(data.*data));
                  % default plotting frame has max variance
  if nargin< 4 | isempty(g.plottimes)
	  g.plottimes = x(plotframes);
  else
	  g.plottimes(find(isnan(g.plottimes))) = x(plotframes);
  end;
  plottimes_set = 1;
end;

if plottimes_set == 1
  ntopos = length(g.plottimes);
  if ntopos > MAX_TOPOS
    fprintf('timtopo(): too many plottimes - only first %d will be shown!\n',MAX_TOPOS);
    g.plottimes = g.plottimes(1:MAX_TOPOS);
    ntopos = MAX_TOPOS;
  end

  if max(g.plottimes) > xmax | min(g.plottimes)< xmin
    fprintf(...
'timtopo(): at least one plottimes value outside of epoch latency range - cannot plot.\n');
    return
  end

  g.plottimes = sort(g.plottimes); % put map latencies in ascending order,
                               % else map lines would cross.
  xshift = [x(2:frames) xmax+1]; % 5/22/2014 Ramon: '+1' was added to avoid errors when time== max(x) in line 231
  plotframes = ones(size(g.plottimes));
  for t = 1:ntopos
    time = g.plottimes(t);
    plotframes(t) = find(time>=x & time < xshift);
  end
end

vlen = length(g.voffsets); % extend voffsets if necessary
i=1;
while vlen< ntopos
        g.voffsets = [g.voffsets g.voffsets(i)];
        i=i+1;
        vlen=vlen+1;
end

%
%%%%%%%%%%%%%%%%  Compute title and axes font sizes %%%%%%%%%%%%%%%
%
pos = get(gca,'Position');
axis('off')
cla % clear the current axes
if pos(4)>0.70
   titlefont= 16;
   axfont = 16;
elseif pos(4)>0.40
   titlefont= 14;
   axfont = 14;
elseif pos(4)>0.30
   titlefont= 12;
   axfont = 12;
elseif pos(4)>0.22
   titlefont= 10;
   axfont = 10;
else
   titlefont= 8;
   axfont = 8;
end

%
%%%%%%%%%%%%%%%% Compute topoplot head width and separation %%%%%%%%%%%%%%%
%
head_sep = 0.2;
topowidth = pos(3)/((6*ntopos-1)/5); % width of each topoplot
if topowidth> 0.25*pos(4) % dont make too large (more than 1/4 of axes width)!
  topowidth = 0.25*pos(4);
end

halfn = floor(ntopos/2);
if rem(ntopos,2) == 1  % odd number of topos
   topoleft = pos(3)/2 - (ntopos/2+halfn*head_sep)*topowidth;
else % even number of topos
   topoleft = pos(3)/2 - ((halfn)+(halfn-1)*head_sep)*topowidth;
end
topoleft = topoleft - 0.01; % adjust left a bit for colorbar

if max(plotframes) > frames |  min(plotframes) < 1
    fprintf('Requested map frame %d is outside data range (1-%d)\n',max(plotframes),frames);
    return
end

%
%%%%%%%%%%%%%%%%%%%% Print times and frames %%%%%%%%%%%%%%%%%%%%%%%%%%
%

fprintf('Scalp maps will show latencies: ');
for t=1:ntopos
  fprintf('%4.0f ',g.plottimes(t));
end
fprintf('\n');
fprintf('                     at frames: ');
for t=1:ntopos
  fprintf('%4d ',plotframes(t));
end
fprintf('\n');

%
%%%%%%%%%%%%%%%%%%%%%%% Plot the data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%% site the plot at bottom of the figure %%%%%%%%%%%%%%%%%%
%
axdata = axes('Units','Normalized','Position',[pos(1) pos(2) pos(3) 0.6*pos(4)],'FontSize',axfont);
set(axdata,'Color',BACKCOLOR);

g.limits = get(axdata,'Ylim');
set(axdata,'GridLineStyle',':')
set(axdata,'Xgrid','off')
set(axdata,'Ygrid','on')
axes(axdata)
axcolor = get(gcf,'Color');
set(axdata,'Color',BACKCOLOR);
pl=plot(x,data(g.plotchans,:)');    % plot the data
if length(g.plotchans)==1
  set(pl,'color','k');
  set(pl,'linewidth',2);
end
xl= xlabel('Latency (ms)');
set(xl,'FontSize',axfont);
yl=ylabel('Potential (\muV)');
set(yl,'FontSize',axfont,'FontAngle','normal');
axis([xmin xmax ymin ymax]);
hold on

%
%%%%%%%%%%%%%%%%%%%%%%%%% Plot zero time line %%%%%%%%%%%%%%%%%%%%%%%%%%
%

if xmin<0 & xmax>0
   plot([0 0],[ymin ymax],'k:','linewidth',1.5);
else
  fprintf('xmin %g and xmax %g do not cross time 0.\n',xmin,xmax)
end
%
%%%%%%%%%%%%%%%%%%%%%%%%% Draw vertical lines %%%%%%%%%%%%%%%%%%%%%%%%%%
%
width  = xmax-xmin;
height = ymax-ymin;
lwidth = 1.5;  % increment line thickness

for t=1:ntopos % dfraw vertical lines through the data at topoplot frames
 if length(g.plotchans)>1 | g.voffsets(t)
  l1 = plot([g.plottimes(t) g.plottimes(t)],...
       [min(data(g.plotchans,plotframes(t))) ...
       g.voffsets(t) + max(data(g.plotchans,plotframes(t)))],'w'); % white underline behind
  set(l1,'linewidth',2);
  l1 = plot([g.plottimes(t) g.plottimes(t)],...
       [min(data(g.plotchans,plotframes(t))) ...
       g.voffsets(t) + max(data(g.plotchans,plotframes(t)))],'b'); % blue line
  set(l1,'linewidth',lwidth);
 end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%% Draw oblique lines %%%%%%%%%%%%%%%%%%%%%%%%%%
%
axall = axes('Position',pos,...
             'Visible','Off','FontSize',axfont);   % whole-gca invisible axes
axes(axall)
set(axall,'Color',BACKCOLOR);
axis([0 1 0 1])
  axes(axall)
  axis([0 1 0 1]);
  set(gca,'Visible','off'); % make whole-figure axes invisible

for t=1:ntopos % draw oblique lines through to the topoplots
  maxdata = max(data(:,plotframes(t))); % max data value at plotframe
  axtp = axes('Units','Normalized','Position',...
       [pos(1)+topoleft+(t-1)*(1+head_sep)*topowidth ...
              pos(2)+0.66*pos(4) ...
                  topowidth ...
                       topowidth*(1+head_sep)]); % this will be the topoplot axes
                       % topowidth]); % this will be the topoplot axes
  axis([-1 1 -1 1]);

  from = changeunits([g.plottimes(t),maxdata],axdata,axall); % data axes
  to   = changeunits([0,-0.74],axtp,axall);                % topoplot axes
  delete(axtp);
  axes(axall);                                             % whole figure axes
  l1 = plot([from(1) to(1)],[from(2) to(2)]);
  set(l1,'linewidth',lwidth);

  hold on
  set(axall,'Visible','off');
  axis([0 1 0 1]);
end
%
%%%%%%%%%%%%%%%%%%%%%%%%% Plot the topoplots %%%%%%%%%%%%%%%%%%%%%%%%%%
%
topoaxes = zeros(1,ntopos);
for t=1:ntopos
       % [pos(3)*topoleft+pos(1)+(t-1)*(1+head_sep)*topowidth ...
  axtp = axes('Units','Normalized','Position',...
       [topoleft+pos(1)+(t-1)*(1+head_sep)*topowidth ...
              pos(2)+0.66*pos(4) ...
                  topowidth topowidth*(1+head_sep)]);
  axes(axtp)                             % topoplot axes
  topoaxes(t) = axtp; % save axes handles
  cla

  if topowidth<0.12
      topoargs = { topoargs{:} 'electrodes' 'off' };
  end
  topoplot( data(:,plotframes(t)),chan_locs, topoargs{:});

  % Else make a 3-D headplot
  %
  % headplot(data(:,plotframes(t)),'chan.spline');

  timetext = [num2str(g.plottimes(t),'%4.0f')];
  % timetext = [num2str(plottimes(t),'%4.0f') ' ms']; % add ' ms'
  text(0.00,0.80,timetext,'FontSize',axfont-3,'HorizontalAlignment','Center'); % ,'fontweight','bold');
end

%
%%%%%%%%%%%%%%%%%%% Plot a topoplot colorbar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
axcb = axes('Position',[pos(1)+pos(3)*0.995 pos(2)+0.62*pos(4) pos(3)*0.02 pos(4)*0.09]);
h=cbar(axcb);                        % colorbar axes
pos_cb = get(axcb,'Position');
set(h,'Ytick',[]);

axes(axall)
set(axall,'Color',axcolor);
%
%%%%%%%%%%%%%%%%%%%%% Plot the color bar '+' and '-' %%%%%%%%%%%%%%%%%%%%%%%%%%
%
text(0.986,0.695,'+','FontSize',axfont,'HorizontalAlignment','Center');
text(0.986,0.625,'-','FontSize',axfont,'HorizontalAlignment','Center');

%
%%%%%%%%%%%%%%%%%%%%%%%%% Plot the plot title if any %%%%%%%%%%%%%%%%%%%%%%%%%%
%
% plot title between data panel and topoplots (to avoid crowding at top of
% figure), on the left
ttl = text(0.03,0.635,g.title,'FontSize',titlefont,'HorizontalAlignment','left'); % 'FontWeight','Bold');

% textent = get(ttl,'extent');
% titlwidth = textent(3);
% ttlpos = get(ttl,'position');
% set(ttl,'position',[     ttlpos(2), ttlpos(3)]);

axes(axall)
set(axall,'layer','top'); % bring component lines to top
for t = 1:ntopos
  set(topoaxes(t),'layer','top'); % bring topoplots to very top
end

  if ~isempty(varargin)
    try,
		if ~isempty( strmatch( 'absmax', varargin))
			text(0.86,0.624,'0','FontSize',axfont,'HorizontalAlignment','Center');
		end;
	catch, end;
  end

%
% Turn on axcopy()
%

% clicking on ERP pop_up topoplot
% -------------------------------
disp('Click on ERP waveform to show scalp map at specific latency');

dat.times = x;
dat.erp   = data;
dat.chanlocs = chan_locs;
dat.options  = topoargs;
dat.srate    = (size(data,2)-1)/(x(end)-x(1))*1000;
dat.axes     = axtp;
dat.line     = l1;

cb_code = [ 'tmppos = get(gca, ''currentpoint'');' ...
            'dattmp = get(gcf, ''userdata'');' ...
            'set(dattmp.line, ''visible'', ''off'');' ...
            'axes(dattmp.axes); cla;' ...
            'latpoint = round((tmppos(1)-dattmp.times(1))/1000*dattmp.srate);' ...
            'topoplot(dattmp.erp(:,latpoint), dattmp.chanlocs, dattmp.options{:});' ...
            'title(sprintf(''%.0f ms'', tmppos(1)));' ...
            'clear latpoint dattmp;' ...
          ];
axcopy;

set(gcf, 'userdata', dat);
set(axdata, 'ButtonDownFcn', cb_code); %windowbuttondownfcn', cb_code);
set(pl, 'ButtonDownFcn', cb_code);

%axcopy(gcf, cb_code);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [gx,dgdx,dgdp] = g_rigid2D(x,P,u,in)
% apply rigid body 2D spatial transform (rotation and translation)
% function [gx,dgdx,dgdp] = g_rigid2D(x,P,u,in)
% IN:
%   - x: [useless]
%   - P: the 3 rigid body transform parameters
%   - u: [useless]
%   - in: contains the original image (in.Y)
% OUT:
%   - gx: vectorized transformed image
%   - dgdx: []
%   - dgdp: gradient of the image wrt to transform parameters

gx = applyRT(P,in);

% the numerical derivation is done here because of interpolation errors:
dgdx = [];
dgdp = zeros(size(gx,1),3);
I = eye(3);
for i=1:3
    dgdp(:,i) = applyRT(P+I(:,i),in) - gx;
end
dgdp = dgdp';

function gx = applyRT(P,in)

tx = P(1);
ty = P(2);
th = P(3);
T = [ 1 0 tx
      0 1 ty
      0 0 1  ];
R = [ cos(th) -sin(th) 0
      sin(th) cos(th)  0
      0       0        1  ];
[nx,ny] = size(in.Y);
[xg,yg] = meshgrid(-nx/2+1:nx/2,-ny/2+1:ny/2);
Txy = T*[xg(:)';yg(:)';ones(1,nx*ny)];
RTxy = R*[Txy(1,:);Txy(2,:);ones(1,nx*ny)];
RTx = reshape(RTxy(1,:)',nx,ny);
RTy = reshape(RTxy(2,:)',nx,ny);
gx = griddata(RTx,RTy,in.Y,xg,yg,'cubic');
gx = gx(:);
gx(isnan(gx)) = 0;

author: MBB-team

repo: VBA-toolbox


text:
function D = prtDistanceHamming(x,y)

[x,y] = prtUtilDistanceParseInputs(x,y);

x = logical(x);
y = logical(y);

%%
%tic
%D2 = sum(bsxfun(@xor,reshape(x,[size(x,1),1,size(x,2)]),reshape(y,[1,size(y,1),size(y,2)])),3);
%toc
%%
D = zeros(size(x,1),size(y,1));
%tic
for iY = 1:size(y,1)
    D(:,iY) = sum(bsxfun(@xor,x,y(iY,:)),2);
end
%toc

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT