text:
function out = isfinite(f)
%ISFINITE   Test if a CLASSICFUN is bounded.
%   ISFINITE(F) returns FALSE if F has any infinite values and TRUE otherwise.
%
% See also ISINF, ISNAN.

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

% Check if the ONEFUN of f is finite:
out = isfinite(f.onefun);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [pitch, yaw] = ueci2angles(reci, veci, ueci)

% convect eci unit vector to rtn angles

% input

%  reci = eci position vector (kilometers)
%  veci = eci velocity vector (kilometers/second)
%  ueci = eci unit vector

% output

%  pitch = pitch angle (radians)
%          positive above the local horizon
%  yaw   = yaw angle (radians)
%          positive in the direction of the angular momentum vector

% Orbital Mechanics with MATLAB

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

% compute radial frame unit vectors

rmag = norm(reci);

xrdl = reci / rmag;

zrdl = cross(reci, veci);

hmag = norm(zrdl);

zrdl = zrdl / hmag;

yrdl = cross(zrdl, xrdl);

% unit vector in radial-tangential-normal frame

umee(1) = dot(ueci, xrdl);

umee(2) = dot(ueci, yrdl);

umee(3) = dot(ueci, zrdl);

% pitch angle (radians)

pitch = asin(umee(1));

% yaw angle (radians)

yaw = atan2(umee(3), umee(2));

author: Sable

repo: mcbench-benchmarks


text:
% StackExchange Signal Processing Q60197
% https://dsp.stackexchange.com/questions/60197
% Implementation of Block Orthogonal Matching Pursuit (BOMP) Algorithm
% References:
%   1.  A
% Remarks:
%   1.  B
% TODO:
% 	1.  C
% Release Notes
% - 1.0.000     19/08/2019
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

numRows         = 8;
% Pay attention: 'numBlocks * numElemntsBlock' should be larger than
% 'numRows'.
numBlocks       = 5;
numElemntsBlock = 3;

paramK = 3;
tolVal = 1e-6;

%% Generate Data

numCols = numBlocks * numElemntsBlock;

mA = randn(numRows, numCols);
vB = randn(numRows, 1);

%% Analysis

% Solving a Problem
vX = SolveLsL0Bomp(mA, vB, numBlocks, paramK, tolVal);

normErr = norm(mA * vX - vB);

disp([' ']);
disp(['Norm of the Error for A x - b - ', num2str(normErr)]);
disp([' ']);

% Comparing to OMP with the number of blocks matches the number of columns
mseErr = mean((SolveLsL0Bomp(mA, vB, numCols, paramK, tolVal) - SolveLsL0Omp(mA, vB, paramK, tolVal)) .^ 2);

disp([' ']);
disp(['MSE of the solution comparison between vanilla OMP to BOMP - ', num2str(mseErr)]);
disp([' ']);

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function [ Acc ] = mi_performance_off2on(CNT,general_initparam)
%CNT{1} = train
%CNT{2} = test
opt = opt_cellToStruct(general_initparam);

% Pre-processing
for onoff=1:length(opt.task)
    CNTch = prep_selectChannels(CNT{onoff}, {'Index', opt.channel_index});
    CNTchfilt =prep_filter(CNTch , {'frequency', opt.band});
    all_SMT{onoff} = prep_segmentation(CNTchfilt, {'interval', opt.time_interval});
    clear CNTch CNTchfilt
end
train = all_SMT{1}; test = all_SMT{2};

% Feature extracion and Classification
[SMT, CSP_W, CSP_D]=func_csp(train,{'nPatterns', opt.CSPFilter});
FT=func_featureExtraction(SMT, {'feature','logvar'});
[CF_PARAM]=func_train(FT,{'classifier','LDA'});

SMT_te=func_projection(test, CSP_W);
FT_te=func_featureExtraction(SMT_te, {'feature','logvar'});
[cf_out]=func_predict(FT_te, CF_PARAM);
[loss out]=eval_calLoss(FT_te.y_dec, cf_out);
Acc=1-loss;
end

author: PatternRecognition

repo: OpenBMI


text:
function varargout = scope_gui(varargin)
% SCOPE_GUI M-file for scope_gui.fig
%      SCOPE_GUI, by itself, creates a new SCOPE_GUI or raises the existing
%      singleton*.
%
%      H = SCOPE_GUI returns the handle to a new SCOPE_GUI or the handle to
%      the existing singleton*.
%
%      SCOPE_GUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in SCOPE_GUI.M with the given input arguments.
%
%      SCOPE_GUI('Property','Value',...) creates a new SCOPE_GUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before scope_gui_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to scope_gui_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help scope_gui

% Last Modified by GUIDE v2.5 24-Sep-2008 14:10:14

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @scope_gui_OpeningFcn, ...
                   'gui_OutputFcn',  @scope_gui_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 scope_gui is made visible.
function scope_gui_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to scope_gui (see VARARGIN)

% Choose default command line output for scope_gui
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes scope_gui wait for user response (see UIRESUME)
% uiwait(handles.figure1);

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

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

% --- Executes on button press in autoscaleButton.
function autoscaleButton_Callback(hObject, eventdata, handles)
% hObject    handle to autoscaleButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

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

% Hints: get(hObject,'String') returns contents of timeBaseEdit as text
%        str2double(get(hObject,'String')) returns contents of timeBaseEdit as a double

% --- Executes during object creation, after setting all properties.
function timeBaseEdit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to timeBaseEdit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

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

% Hints: get(hObject,'String') returns contents of vscaleEdit as text
%        str2double(get(hObject,'String')) returns contents of vscaleEdit as a double

% --- Executes during object creation, after setting all properties.
function vscaleEdit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to vscaleEdit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on slider movement.
function timeBaseSlider_Callback(hObject, eventdata, handles)
% hObject    handle to timeBaseSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

% --- Executes during object creation, after setting all properties.
function timeBaseSlider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to timeBaseSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function vscaleSlider_Callback(hObject, eventdata, handles)
% hObject    handle to vscaleSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

% --- Executes during object creation, after setting all properties.
function vscaleSlider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to vscaleSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

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

% Hints: get(hObject,'String') returns contents of updateIntervalEdit as text
%        str2double(get(hObject,'String')) returns contents of updateIntervalEdit as a double

% --- Executes during object creation, after setting all properties.
function updateIntervalEdit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to updateIntervalEdit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


text:
function varargout = kj_grid(varargin)
% KJ_GRID Application M-file for kj_grid.fig
%    FIG = KJ_GRID launch kj_grid GUI.
%    KJ_GRID('callback_name', ...) invoke the named callback.

% Last Modified by GUIDE v2.0 12-Jun-2002 11:51:21

if nargin == 1  % LAUNCH GUI

  fig = openfig(mfilename,'reuse');

  % Use system color scheme for figure:
  set(fig,'Color',get(groot,'defaultUicontrolBackgroundColor'));

  % Generate a structure of handles to pass to callbacks, and store it.
  handles = guihandles(fig);
  guidata(fig, handles);

  if varargin{1} == 1
    set(handles.lblSpacingDepth, 'String', 'Latitude spacing [deg]:');
    set(handles.txtSpacingDepth, 'String', '0.1');
    set(handles.lblSpacingHorizontal, 'String', 'Longitude spacing [deg]:');
    set(handles.txtSpacingHorizontal, 'String', '0.1');
  else
    set(handles.lblSpacingDepth, 'String', 'Depth spacing [km]:');
    set(handles.txtSpacingDepth, 'String', '1');
    set(handles.lblSpacingHorizontal, 'String', 'Horizontal spacing [km]:');
    set(handles.txtSpacingHorizontal, 'String', '1');
  end

  if get(handles.chkRandomNode, 'Value') == 1
    set(handles.txtNumberCalculationNode, 'Enable', 'on');
    set(handles.chkRandomArea, 'Enable', 'on');
    if get(handles.chkRandomArea, 'Value') == 1
      set(handles.txtNumberCalculationArea, 'Enable', 'on');
    else
      set(handles.txtNumberCalculationArea, 'Enable', 'off');
    end
  else
    set(handles.txtNumberCalculationNode, 'Enable', 'off');
    set(handles.chkRandomArea, 'Enable', 'off');
    set(handles.txtNumberCalculationArea, 'Enable', 'off');
  end

  if get(handles.chkLearningPeriod, 'Value') == 1
    set(handles.txtLearningPeriod, 'Enable', 'on');
  else
    set(handles.txtLearningPeriod, 'Enable', 'off');
  end

  if get(handles.chkForecastPeriod, 'Value') == 1
    set(handles.txtForecastPeriod, 'Enable', 'on');
  else
    set(handles.txtForecastPeriod, 'Enable', 'off');
  end

  if get(handles.radNumber, 'Value') == 1
    set(handles.txtNumber, 'Enable', 'on');
    set(handles.radRadius, 'Value', [0]);
    set(handles.txtRadius, 'Enable', 'off');
  else
    set(handles.txtNumber, 'Enable', 'off');
    set(handles.radRadius, 'Value', [1]);
    set(handles.txtRadius, 'Enable', 'on');
  end

  if get(handles.chkSaveParameter, 'Value') == 1
    set(handles.btnSaveParameter, 'Enable', 'on');
  else
    set(handles.btnSaveParameter, 'Enable', 'off');
  end

  if get(handles.chkMinMag, 'Value') == 1
    set(handles.txtMinMag, 'Enable', 'off');
  else
    set(handles.txtMinMag, 'Enable', 'on');
  end

  % Wait for callbacks to run and window to be dismissed:
  uiwait(fig);

  if nargout > 0
    varargout{1} = fig;
  end

elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

  try
    [varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
  catch
    disp(lasterr);
  end

end

% --------------------------------------------------------------------
function varargout = radNumber_Callback(h, eventdata, handles, varargin)

if get(handles.radNumber, 'Value') == 1
  set(handles.txtNumber, 'Enable', 'on');
  set(handles.radRadius, 'Value', [0]);
  set(handles.txtRadius, 'Enable', 'off');
  set(handles.radRectangle, 'Value', [0]);
  set(handles.txtSizeRectHorizontal, 'Enable', 'off');
  set(handles.txtSizeRectDepth, 'Enable', 'off');
else
  set(handles.txtNumber, 'Enable', 'off');
  set(handles.radRadius, 'Value', [1]);
  set(handles.txtRadius, 'Enable', 'on');
end

% --------------------------------------------------------------------
function varargout = radRadius_Callback(h, eventdata, handles, varargin)

if get(handles.radRadius, 'Value') == 1
  set(handles.txtRadius, 'Enable', 'on');
  set(handles.radNumber, 'Value', [0]);
  set(handles.txtNumber, 'Enable', 'off');
  set(handles.radRectangle, 'Value', [0]);
  set(handles.txtSizeRectHorizontal, 'Enable', 'off');
  set(handles.txtSizeRectDepth, 'Enable', 'off');
else
  set(handles.txtRadius, 'Enable', 'off');
  set(handles.radNumber, 'Value', [1]);
  set(handles.txtNumber, 'Enable', 'on');
end

% --------------------------------------------------------------------
function varargout = radRectangle_Callback(h, eventdata, handles, varargin)

if get(handles.radRectangle, 'Value') == 1
  set(handles.txtSizeRectHorizontal, 'Enable', 'on');
  set(handles.txtSizeRectDepth, 'Enable', 'on');
  set(handles.radRadius, 'Value', [0]);
  set(handles.txtRadius, 'Enable', 'off');
  set(handles.radNumber, 'Value', [0]);
  set(handles.txtNumber, 'Enable', 'off');
else
  set(handles.txtRadius, 'Enable', 'off');
  set(handles.radNumber, 'Value', [1]);
  set(handles.txtNumber, 'Enable', 'on');
end

% --------------------------------------------------------------------
function varargout = chkRandomNode_Callback(h, eventdata, handles, varargin)

if get(handles.chkRandomNode, 'Value') == 1
  set(handles.txtNumberCalculationNode, 'Enable', 'on');
  set(handles.chkRandomArea, 'Enable', 'on');
  if get(handles.chkRandomArea, 'Value') == 1
    set(handles.txtNumberCalculationArea, 'Enable', 'on');
  else
    set(handles.txtNumberCalculationArea, 'Enable', 'off');
  end
else
  set(handles.txtNumberCalculationNode, 'Enable', 'off');
  set(handles.chkRandomArea, 'Enable', 'off');
  set(handles.txtNumberCalculationArea, 'Enable', 'off');
end

% --------------------------------------------------------------------
function varargout = chkRandomArea_Callback(h, eventdata, handles, varargin)

if get(handles.chkRandomArea, 'Value') == 1
  set(handles.txtNumberCalculationArea, 'Enable', 'on');
else
  set(handles.txtNumberCalculationArea, 'Enable', 'off');
end

% --------------------------------------------------------------------
function varargout = chkLearningPeriod_Callback(h, eventdata, handles, varargin)

if get(handles.chkLearningPeriod, 'Value') == 1
  set(handles.txtLearningPeriod, 'Enable', 'on');
else
  set(handles.txtLearningPeriod, 'Enable', 'off');
end

% --------------------------------------------------------------------
function varargout = chkForecastPeriod_Callback(h, eventdata, handles, varargin)

if get(handles.chkForecastPeriod, 'Value') == 1
  set(handles.txtForecastPeriod, 'Enable', 'on');
else
  set(handles.txtForecastPeriod, 'Enable', 'off');
end

% --------------------------------------------------------------------
function varargout = chkSaveParameter_Callback(h, eventdata, handles, varargin)

if get(handles.chkSaveParameter, 'Value') == 1
  set(handles.btnSaveParameter, 'Enable', 'on');
else
  set(handles.btnSaveParameter, 'Enable', 'off');
end

% --------------------------------------------------------------------
function varargout = chkMinMag_Callback(h, eventdata, handles, varargin)

if get(handles.chkMinMag, 'Value') == 1
  set(handles.txtMinMag, 'Enable', 'off');
else
  set(handles.txtMinMag, 'Enable', 'on');
end

% --------------------------------------------------------------------
function varargout = btnSaveParameter_Callback(h, eventdata, handles, varargin)

[newfile, newpath] = uiputfile('*.mat','Choose parameter file');
% Cancel pressed?
if isequal(newfile, 0)  ||  isequal(newpath, 0)
  return;
end
% Everything ok?
newfile = [newpath newfile];
if length(newfile) > 1
  set(handles.lblSaveParameter, 'String', newfile);
end

% --------------------------------------------------------------------
function varargout = btnCancel_Callback(h, eventdata, handles, varargin)

handles.answer = 0;
guidata(h, handles);
uiresume(handles.figGrid);

% --------------------------------------------------------------------
function varargout = btnOk_Callback(h, eventdata, handles, varargin)

% Validation of inputs
if get(handles.radNumber, 'Value') == 1
  Tmp = str2double(get(handles.txtNumber, 'String'));
  if (isnan(Tmp) | (Tmp ~= round(Tmp)) | (Tmp <= 0))
    errordlg('Number of events must be a positive integer value');
    return;
  end
end

if get(handles.radRadius, 'Value') == 1
  Tmp = str2double(get(handles.txtRadius, 'String'));
  if (isnan(Tmp) | (Tmp <= 0))
    errordlg('Constant radius must be a positive float value');
    return;
  end
end

Tmp = str2double(get(handles.txtMinimumNumber, 'String'));
if (isnan(Tmp) | (Tmp ~= round(Tmp)) | (Tmp <= 0))
  errordlg('Minimum number of events must be a positive integer value');
  return;
end

Tmp = str2double(get(handles.txtSpacingHorizontal, 'String'));
if (isnan(Tmp) | (Tmp <= 0))
  errordlg('Horizontal spacing must be a positive float value');
  return;
end

Tmp = str2double(get(handles.txtSpacingDepth, 'String'));
if (isnan(Tmp) | (Tmp <= 0))
  errordlg('Depth spacing must be a positive float value');
  return;
end

if get(handles.chkLearningPeriod, 'Value') == 1
  Tmp = str2double(get(handles.txtLearningPeriod, 'String'));
  if (isnan(Tmp) | (Tmp <= 0))
    errordlg('Learning period must be a positive float value');
    return;
  end
end

if get(handles.chkForecastPeriod, 'Value') == 1
  Tmp = str2double(get(handles.txtForecastPeriod, 'String'));
  if (isnan(Tmp) | (Tmp <= 0))
    errordlg('Forecasting period must be a positive float value');
    return;
  end
end

Tmp = str2double(get(handles.txtSplitTime, 'String'));
if (isnan(Tmp) | (Tmp <= 0))
  errordlg('Splitting time must be a positive float value');
  return;
end

if get(handles.chkRandomNode, 'Value') == 1
  Tmp = str2double(get(handles.txtNumberCalculationNode, 'String'));
  if (isnan(Tmp) | (Tmp ~= round(Tmp)) | (Tmp <= 0))
    errordlg('Number of calculations per node must be a positive integer value');
    return;
  end
  if get(handles.chkRandomArea, 'Value') == 1
    Tmp = str2double(get(handles.txtNumberCalculationArea, 'String'));
    if (isnan(Tmp) | (Tmp ~= round(Tmp)) | (Tmp <= 0))
      errordlg('Number of calculations over area must be a positive integer value');
      return;
    end
  end
end

if get(handles.chkSaveParameter, 'Value') == 1
  Tmp = get(handles.lblSaveParameter, 'String');
  if isempty(Tmp)
    errordlg('Please select a filename for storing the parameters.');
    return;
  end
end

if get(handles.chkMinMag, 'Value') == 0
  Tmp = str2double(get(handles.txtMinMag, 'String'));
  if (isnan(Tmp) | ((Tmp * 10) ~= round(Tmp * 10)) |(Tmp <= 0))
    errordlg('Minimum magnitude must be a positive float value (precision 0.1)');
    return;
  end
end

Tmp = str2double(get(handles.txtMaxMag, 'String'));
if (isnan(Tmp) | ((Tmp * 10) ~= round(Tmp * 10)) |(Tmp <= 0))
  errordlg('Maximum magnitude must be a positive float value (precision 0.1)');
  return;
end

% OK, proceed
handles.answer = 1;
guidata(h, handles);
uiresume(handles.figGrid);

author: CelsoReyes


text:
function test_suite=test_cross_neighborhood()
% tests for cosmo_cross_neighborhood
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_cross_neighborhood_time_freq()
    helper_test_cross_neighborhood(false);

function test_cross_neighborhood_chan_time_freq()
    if cosmo_skip_test_if_no_external('fieldtrip')
        return;
    end
    helper_test_cross_neighborhood(true);

function helper_test_cross_neighborhood(can_use_chan_nbrhood)
    ds_full=cosmo_synthetic_dataset('type','timefreq','size','big');
    msk=cosmo_match(ds_full.fa.chan,@(x)x<20);
    ds_full=cosmo_slice(ds_full,msk,2);
    ds_full=cosmo_dim_prune(ds_full);

    % make it sparse
    nfeatures=size(ds_full.samples,2);
    nkeep=round(nfeatures/4);
    fdim_values=ds_full.a.fdim.values;
    nchan=numel(fdim_values{1});
    nfreq=numel(fdim_values{2});
    ntime=numel(fdim_values{3});

    while true
        rp=randperm(nfeatures);
        ds=cosmo_slice(ds_full,rp(1:nkeep),2);
        ds=cosmo_dim_prune(ds);
        n=numel(ds.a.fdim.values{1});
        ds.a.fdim.values{1}=ds.a.fdim.values{1}(randperm(n));

        if numel(unique(ds.fa.chan))==nchan && ...
                    numel(unique(ds.fa.freq))==nfreq && ...
                    numel(unique(ds.fa.time))==ntime
            break
        end
    end

    nfeatures=size(ds.samples,2);

    % define neighborhoods
    freq_nbrhood=cosmo_interval_neighborhood(ds,'freq','radius',2);
    time_nbrhood=cosmo_interval_neighborhood(ds,'time','radius',1);
    if can_use_chan_nbrhood
        chan_nbrhood=cosmo_meeg_chan_neighborhood(ds,'count',5,...
                                'chantype','all','label','dataset');
    else
        chan_nbrhood='dummy';
    end

    all_nbrhoods={chan_nbrhood, freq_nbrhood, time_nbrhood};
    ndim=numel(all_nbrhoods);
    dim_labels={'chan';'freq';'time'};

    ntest=5;  % number of positions to test

    for i=7:-1:1
        use_chan=i<=4;
        use_freq=mod(i,2)==1;
        use_time=mod(ceil(i/2),2)==1;

        if ~can_use_chan_nbrhood && use_chan
            % no support for channel neighborhood, skip
            continue;
        end

        use_dim_msk=[use_chan;use_freq;use_time];
        nbrhood=cosmo_cross_neighborhood(ds,all_nbrhoods(use_dim_msk),...
                                                'progress',false);
        assertEqual(nbrhood.a.fdim.labels,dim_labels(use_dim_msk));
        assertEqual(fieldnames(nbrhood.fa),dim_labels(use_dim_msk));

        n=numel(nbrhood.neighbors);
        rp=randperm(n);

        for iter=1:min(n,ntest)
            pos=rp(iter);
            % verify neighborhoods in ds

            ds_fa=cosmo_slice(ds.fa,nbrhood.neighbors{pos},2,'struct');
            nbr_fa=cosmo_slice(nbrhood.fa,pos,2,'struct');

            nbr_msk=true(1,nfeatures);

            for dim=1:ndim
                dim_label=dim_labels{dim};
                if use_dim_msk(dim)
                    dim_nbrhood=all_nbrhoods{dim};

                    j=find(dim_nbrhood.fa.(dim_label)==nbr_fa.(dim_label));
                    assert(numel(j)==1);

                    m=false(1,nfeatures);
                    m(dim_nbrhood.neighbors{j})=true;
                else
                    assert(~isfield(nbrhood.fa,dim_label))
                    m=true(1,nfeatures);
                end

                nbr_msk=nbr_msk & m;

                fa=cosmo_slice(ds.fa.(dim_label),m,2);
                assert(isempty(setdiff(ds_fa.(dim_label),fa)));
            end

            assertEqual(nbrhood.neighbors{pos},find(nbr_msk));

            % test agreement between the crossed nbrhood and the
            % individual neighborhoods
            dim_nbr_msk=true(1,nfeatures);
            dim_pos=0;
            for dim=1:ndim
                if ~use_dim_msk(dim)
                    continue;
                end
                dim_pos=dim_pos+1;
                dim_label=dim_labels{dim};
                nbr_fa=nbrhood.fa.(dim_label);
                nbr_values=nbrhood.a.fdim.values{dim_pos}(nbr_fa(pos));

                dim_nbrhood=all_nbrhoods{dim};
                dim_nbr_values=dim_nbrhood.a.fdim.values{1};
                nbr_msk=cosmo_match(dim_nbr_values,nbr_values);

                m=false(size(dim_nbr_msk));
                m(dim_nbrhood.neighbors{nbr_msk})=true;
                dim_nbr_msk=dim_nbr_msk & m;
            end
            assertEqual(nbrhood.neighbors{pos},find(dim_nbr_msk));
        end
    end

function test_cross_neighborhood_transpose
    opt=struct();
    opt.progress=false;
    ds=cosmo_synthetic_dataset('type','timefreq','size','normal');
    ds=cosmo_dim_remove(ds,'chan');

    nh_time=cosmo_interval_neighborhood(ds,'time','radius',1);
    nh_freq=cosmo_interval_neighborhood(ds,'freq','radius',0);

    nh=cosmo_cross_neighborhood(ds,{nh_freq,nh_time},opt);

    cp=cosmo_cartprod(repmat({[false,true]},4,1));
    n=size(cp,1);

    for k=1:n
        t_label=cp{k,1};
        t_value=cp{k,2};
        t_elem1=cp{k,3};
        t_elem2=cp{k,4};

        ds2=ds;

        if t_label
            ds2.a.fdim.labels=ds2.a.fdim.labels';
        end

        if t_value
            ds2.a.fdim.values=ds2.a.fdim.values';
        end

        if t_elem1
            ds2.a.fdim.values{1}=ds2.a.fdim.values{1}';
        end

        if t_elem2
            ds2.a.fdim.values{2}=ds2.a.fdim.values{2}';
        end

        nh2_time=cosmo_interval_neighborhood(ds2,'time','radius',1);
        nh2_freq=cosmo_interval_neighborhood(ds2,'freq','radius',0);

        nh2=cosmo_cross_neighborhood(ds2,{nh2_freq,nh2_time},opt);
        assertEqual(nh2.a,nh.a);
        assertEqual(nh2.fa,nh.fa);
        assertEqual(nh2.neighbors,nh.neighbors);
    end

function test_cross_neighborhood_exceptions()
    ds=cosmo_synthetic_dataset('type','meeg','size','big');
    time_nbrhood=cosmo_interval_neighborhood(ds,'time','radius',1);

    % test exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                        cosmo_cross_neighborhood(varargin{:}),'');

    aet(ds,{});
    aet(ds,struct);
    aet(ds,struct,time_nbrhood);
    aet(ds,time_nbrhood,time_nbrhood);
    aet(ds,ds);

    % should be fine
    assertEqual(cosmo_cross_neighborhood(ds,{time_nbrhood}),time_nbrhood);

    % no values too big
    time_nbrhood2=time_nbrhood;
    time_nbrhood2.neighbors{1}=1e9;
    aet(ds,{time_nbrhood2});

    % non-integers not supported
    time_nbrhood2=time_nbrhood;
    time_nbrhood2.neighbors{1}=1.5;
    aet(ds,{time_nbrhood2});

    % illegal labels
    time_nbrhood2=time_nbrhood;
    time_nbrhood2.a.fdim.labels{1}='foo';
    time_nbrhood2.fa.foo=time_nbrhood.fa.time;
    aet(ds,{time_nbrhood2});

    % duplicate labels
    aet(ds,{time_nbrhood,time_nbrhood});

function test_cross_neighborhood_unsorted_neighbors
    ds=cosmo_synthetic_dataset();
    nh=cosmo_interval_neighborhood(ds,'i','radius',1);

    nh_unsorted=nh;
    nh_unsorted.neighbors=cellfun(@(x)x(randperm(numel(x))),...
                                nh_unsorted.neighbors,...
                                'UniformOutput',false);
    assertEqual(nh,cosmo_cross_neighborhood(ds,{nh_unsorted}));

function test_cross_neighborhood_progress()
    if cosmo_skip_test_if_no_external('!evalc')
        return;
    end

    ds=cosmo_synthetic_dataset();
    nh1=cosmo_interval_neighborhood(ds,'i','radius',0);
    nh2=cosmo_interval_neighborhood(ds,'j','radius',0);
    f=@()cosmo_cross_neighborhood(ds,{nh1,nh2});
    res=evalc('f();');
    assert(~isempty(strfind(res,'[####################]')));
    assert(~isempty(strfind(res,'crossing neighborhoods')));

function test_warning_weird_dimension_order
    if cosmo_skip_test_if_no_external('fieldtrip')
        return;
    end

    ds=cosmo_synthetic_dataset('type','timefreq','size','big');
    keep=ds.fa.chan<20;
    ds=cosmo_slice(ds,keep,2);
    ds=cosmo_dim_prune(ds);

    nh_time=cosmo_interval_neighborhood(ds,'time','radius',0);
    nh_freq=cosmo_interval_neighborhood(ds,'freq','radius',0);
    nh_chan=cosmo_meeg_chan_neighborhood(ds,'count',4,...
                                    'chantype','meg_planar',...
                                    'label','dataset');

    nh_cell={nh_chan,nh_freq,nh_time};
    for ndim=1:3
        idxs=nchoosek(1:3,ndim);
        order=perms(1:ndim);

        for k=1:size(idxs,1)
            for j=1:size(order,1)

                idx=idxs(k,order(j,:));
                nh_sel=nh_cell(idx);

                is_weird_order=~issorted(idx);

                func=@()cosmo_cross_neighborhood(ds,nh_sel,...
                                    'progress',false);
                assert_warning_being_shown(func,is_weird_order);
            end
        end
    end

function assert_warning_being_shown(func, flag)
    warning_state=cosmo_warning();
    cleaner=onCleanup(@()cosmo_warning(warning_state));

    cosmo_warning('reset');
    cosmo_warning('off');

    % initially must be empty
    s=cosmo_warning();
    assertEqual(s.shown_warnings,{});

    func();

    s=cosmo_warning();
    if flag
        % must have warning added
        assertTrue(iscellstr(s.shown_warnings));
        assertTrue(numel(s.shown_warnings)>0);
    else
        % still without warning
        assertEqual(s.shown_warnings,{});
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function hcell_node_display ( node_xy_file_name )

%*****************************************************************************80
%
%% HCELL_NODE_DISPLAY displays the nodes in the H-Cell.
%
%  Usage:
%
%    hcell_node_display ( node_xy_file_name )
%
%    A typical invocation might be
%
%      hcell_node_display ( 'xy3.txt' )
%
%    or
%
%      hcell_node_display ( 'xy6.txt' )
%
%    But if you simply say
%
%      hcell_node_display
%
%    the program will give you a chance to enter the file names
%    interactively.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, filename NODE_XY_FILE_NAME, the name of a file
%    containing the XY coordinates of nodes.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HCELL_NODE_DISPLAY:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Display the nodes in the HCELL problem.\n' );
%
%  Do we have the XY file?
%
  if ( nargin < 1 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'HCELL_NODE_DISPLAY:\n' );
    node_xy_file_name = input ( 'Enter the name of the XY coordinate file:' );

  end

  [ m1, n1 ] = r8mat_header_read ( node_xy_file_name );

  node_xy = r8mat_data_read ( node_xy_file_name, m1, n1 );
%
%  Plot the nodes.
%
  plot ( node_xy(1,:), node_xy(2,:), 'bo' )
%
%  Add the boundary of the region to the plot.
%
  hcell_boundary_add ( 'r' );
%
%  Add the invisible bounding box.
%
  hcell_box_add ( 'w' )

  axis equal

  title ( 'Node coordinates' )

  text ( 45.0, 20.0, node_xy_file_name )

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HCELL_NODE_DISPLAY:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
    column_num = -1;
    return;
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return;
  end

  column_num = s_word_count ( line );

  return
end
function row_num = file_row_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_ROW_COUNT counts the number of row records in a file.
%
%  Discussion:
%
%    Each input line is a "RECORD".
%
%    The records are divided into three groups:
%
%    * BLANK LINES (nothing but blanks)
%    * COMMENT LINES (begin with a '#')
%    * DATA RECORDS (anything else)
%
%    The value returned by the function is the number of data records.
%
%    By the way, if the MATLAB routine FGETS is used, instead of
%    FGETL, then the variable LINE will include line termination
%    characters, which means that a blank line would not actually
%    have zero characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer ROW_NUM, the number of rows found.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_ROW_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_ROW_COUNT - Error!' );
    row_num = -1;
    return;
  end

  blank_num = 0;
  comment_num = 0;
  row_num = 0;

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    record_num = record_num + 1;
    record_length = s_len_trim ( line );

    if ( record_length <= 0 )
      blank_num = blank_num + 1;
    elseif ( line(1) == '#' )
      comment_num = comment_num + 1;
    else
      row_num = row_num + 1;
    end

  end

  fclose ( input_unit );

  return
end
function hcell_box_add ( box_color )

%*****************************************************************************80
%
%% HCELL_BOX_ADD adds a box around the HCELL region.
%
%  Discussion:
%
%    This routine can only ADD graphics to an already existing plot.
%    So call PLOT first and then this routine.
%
%    The box is normally drawn with the color white ('w').
%    In this case, it doesn't actually appear on the screen.
%    However, it does force MATLAB to make the picture bigger,
%    giving us space to write some titles.
%
%  Usage:
%
%    hcell_box_add ( box_color )
%
%    A typical invocation might be
%
%      hcell_box_add ( 'w' )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character BOX_COLOR, the color to be used to draw the box.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HCELL_BOX_ADD:\n' );
  fprintf ( 1, '  Add a bounding box to a plot of the HCELL problem.\n' );
%
%  Do we have a color?
%
  if ( nargin < 1 )
    box_color = 'w'
  end
%
%  Set the coordinates of the invisible bounding box that
%  allows us to display the file name within the plot area.
%
  box_x = [ -5.00,  110.00, 110.00, -5.00, -5.00 ];
  box_y = [ -5.00,   -5.00, 21.00,  21.00, -5.00 ];

  line ( box_x, box_y, 'color', box_color )

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HCELL_BOX_ADD:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function table = r8mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%
  table = [];
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %f' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the file.\n' );
    error ( 'R8MAT_DATA_READ - Error!' );
    return;
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = r8mat_header_read ( input_filename )

%*****************************************************************************80
%
%% R8MAT_HEADER_READ reads the header from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function modelPlane = sr_extract_plane(imgPath, imgName, opt)

% SR_EXTRACT_PLANE:
%
% Extract plane model from an image
%
% The code is adapted and modified from the paper
%
% Jia-Bin Huang, Sing Bing Kang, Narendra Ahuja, Johannes Kopf
% Image Completion using Planar Structure Guidance
% ACM Transactions on Graphics (Proceedings of SIGGRAPH 2014).
%
% Output: modelPlane
% The model plane has the following data structure
%
% modelPlane
% modelPlane.vpData:                        Vanishing point data
% modelPlane.numPlane                       Number of planes
% modelPlane.plane{indPlane}.vLine          The vanishing line of the plane
% modelPlane.plane{indPlane}.imgPlaneProb;  The planar location density
% modelPlane.plane{indPlane}.sourceVP       Which two VPs form the plane
% modelPlane.plane{indPlane}.rotPar(vpInd)  Rotation parameters for aligning
%                                           two sets of lines with the x-axis                                           x-axis
% modelPlane.plane{indPlane}.postProb       Posteior probability of the plane

% =========================================================================
% Vanishing point detection
% =========================================================================

vpFilePath = 'cache\vpdetection';
vpFileName  = [imgName(1:end-4), '-vanishingpoints.txt'];

recomputeFlag = 1;
if(~exist(fullfile(vpFilePath, 'text', vpFileName), 'file') || recomputeFlag)
    vpExeFile = 'source\EdgeDetectTest.exe';
    vpDetectCMD = [vpExeFile, ' -indir ', imgPath, ' -infile ', imgName, ' -outdir ', vpFilePath];
    system(vpDetectCMD);
end
% Read vanishing point data
vpData = sr_read_vpdata(fullfile(vpFilePath, 'text', vpFileName));

img = imread(fullfile(imgPath, imgName));

% =========================================================================
% Plane localization
% =========================================================================

modelPlane = sr_detect_plane_from_vp(vpData, img, opt);
% figure(1), imshow(modelPlane.postProb(:,:,1:3));

visVPFlag = 0;
if(visVPFlag)
    [vpVis, planeVis] = vis_vp(img, vpData, modelPlane.postProb(:,:,1:3));
end
end

function [vpVis, planeVis] = vis_vp(img, vpData, postProb)

img = im2double(img);

[imgH, imgW, nCh] = size(img);
vpVis = zeros(imgH, imgW, nCh);
planeVis = zeros(imgH, imgW, nCh);
whiteImg = ones(imgH, imgW, nCh);
% blackImg = zeros()
alphaW = 0.5;

imgW = img*(1 - alphaW) + whiteImg*(alphaW);

% Plot VP
h1 = figure(1);
imshow(imgW); hold on;
lineWidth = 4;
for vpInd = 1: vpData.numVP
    vpCurr = vpData.vp{vpInd};
    for lineInd = 1: vpCurr.numLines
        lineCur = vpCurr.lines(lineInd,:);
        lineCur = lineCur + 1;
        if(vpInd == 1)
            plot(lineCur([1,3]), lineCur([2,4]), 'r', 'LineWidth', lineWidth);
        elseif(vpInd == 2)
            plot(lineCur([1,3]), lineCur([2,4]), 'g', 'LineWidth', lineWidth);
        elseif(vpInd == 3)
            plot(lineCur([1,3]), lineCur([2,4]), 'b', 'LineWidth', lineWidth);
        end
    end
end
hold off;
print(h1, '-dpng', fullfile('paper\planar_struct_SR_CVPR2015\figures\plane', 'vpdetection.png'));

% Plot posterior
alphaP = 0.9;
imgP = (1 - alphaP)*img + alphaP*postProb;
figure(2); imshow(imgP);
imwrite(imgP, fullfile('paper\planar_struct_SR_CVPR2015\figures\plane', 'planedetection.png'));

end
function vpData = sr_read_vpdata(fileName)

% SC_READ_VPDATA: read the data from vanishing point detection algorithm
% Input:
%   - fileName: the txt file containing the pre-computed vanishing point
%   detection code
% Output:
%   - vpData
%   The data structure of vpData
%       - vpData.numVP: number of detected vanishing points
%       - vpData.vp{i}.pos: the vanishing point position in the homogenous coordiante
%       - vpData.vp{i}.score: the score of the vanishing point
%       - vpData.vp{i}.numLines: number of lines supporting the vanishing point
%       - vpData.vp{i}.lines{j}.p1: (x1, y1): starting position
%       - vpData.vp{i}.lines{j}.p2: (x2, y2): ending position
%       - vpData.vp{i}.lines{j}.length: length of the line segment

vpData = [];

% Read data
fid = fopen(fileName);

%% Parse VP positions
temp = fscanf(fid, '%s ', [1 5]);
numVP = 0;
readVPFlag = 1;
VP = [];
while(readVPFlag)
    numVP = numVP + 1;
    vpCurr = fscanf(fid, '%g %g %g %g %g', [5 1]);
    if(~isempty(vpCurr))
        VP(:,numVP) = vpCurr;
    else
        temp = fscanf(fid, '%s ', [1 6]);
        readVPFlag = 0;
    end
end
VP = VP';

vpData.numVP = size(VP, 1);

% Save VP position data
for i = 1: vpData.numVP
    vpData.vp{i}.pos = VP(i, 1:3);
    vpData.vp{i}.score = VP(i, 4);
    vpData.vp{i}.numLines = VP(i, 5);
end

%% Parse each set of line segments for the corresponding VP
for i = 1: vpData.numVP
    numLine = fscanf(fid, '%d ', [1 1]);
    lines = fscanf(fid, '%g %g %g %g %g', [5 numLine]);
    vpData.vp{i}.lines = lines';
end

fclose(fid);
end

function modelPlane = sr_detect_plane_from_vp(vpData, img, opt)

% SC_DETECT_PLANE_FROM_VP: simple plane detection algorithm
% Input:
%     - vpData: vanishing point data
%     - img: input image
%     - mask: hole mask
% Output:
%     - modelPlane

%%

modelPlane = [];

% === Setting up ===
[imgH, imgW, ch] = size(img);
HfilterX = fspecial('gaussian', [1, opt.filterSize], opt.filterSigma);
HfilterY = HfilterX';
% fspecial('gaussian', opt.filterSize, opt.filterSigma);

img = im2double(img);

% === Supporting lines spatial support estimation ===
shapeInserter = vision.ShapeInserter('Shape', 'Lines','BorderColor', 'White');

for i = 1: vpData.numVP
    % The support lines
    imgLines = zeros(imgH, imgW);
    imgLines = step(shapeInserter, imgLines, int16(round(vpData.vp{i}.lines(:,1:4))));
    % Spatial density estimation via blurring
    imgLinesPosMap = imgLines;
    for k = 1:opt.numFilterIter
        imgLinesPosMap = imfilter(imgLinesPosMap, HfilterX, 'conv', 'replicate');
    end
    for k = 1:opt.numFilterIter
        imgLinesPosMap = imfilter(imgLinesPosMap, HfilterY, 'conv', 'replicate');
    end

    % Save results
    modelPlane.vp{i}.imgLines = imgLines;
    modelPlane.vp{i}.imgLinesPosMap = imgLinesPosMap;
end

% === Estimate plane support and plane parameters ===
numPlane = (vpData.numVP)*(vpData.numVP-1)/2;
% Initialize plane data
modelPlane.plane = cell(numPlane, 1);

indPlane = 1;
% A pair of vanishing points forms a plane hypothesis
for i = 1: vpData.numVP - 1
    for j = i+1: vpData.numVP
        % Compute the vanishing line
        modelPlane.plane{indPlane}.vLine = vLineFromTwoVP(vpData.vp{i}.pos, vpData.vp{j}.pos);
        % Element-wise product of two support line density
        modelPlane.plane{indPlane}.imgPlaneProb = modelPlane.vp{i}.imgLinesPosMap.*modelPlane.vp{j}.imgLinesPosMap; % Product of two probability maps

        %         modelPlane.plane{indPlane}.imgPlaneProb(mask) = 1e-10;
        modelPlane.plane{indPlane}.sourceVP = [i, j];

        indPlane = indPlane + 1;
    end
end

% === Compute rectified rotation parameters ===

for i = 1: numPlane
    for vpInd = 1: 2

        linesCurr = vpData.vp{modelPlane.plane{i}.sourceVP(vpInd)}.lines;
        invalidLineInd = linesCurr(:,5) == 0;
        linesCurr = linesCurr(~invalidLineInd,:);
        numLines = size(linesCurr, 1);

        vLineCurr = modelPlane.plane{i}.vLine;

        % Rectified homography
        H = eye(3);
        H(3,:) = vLineCurr;

    end
end

% === Add a fronto-parallel plane ===

modelPlane.plane{indPlane}.vLine = [0 0 1];
modelPlane.plane{indPlane}.imgPlaneProb = opt.fpPlaneProb*ones(imgH, imgW);
modelPlane.plane{indPlane}.score = sum(modelPlane.plane{indPlane}.imgPlaneProb(:));

numPlane = numPlane + 1;

modelPlane.numPlane = numPlane;

% === Compute posterior probability ===

planeProb = zeros(imgH, imgW, numPlane);
for i = 1 : numPlane
    planeProb(:,:,i) = modelPlane.plane{i}.imgPlaneProb;
end
planeProbSum = sum(planeProb, 3);
planeProb = bsxfun(@rdivide, planeProb, planeProbSum);
planeProb = planeProb + 0.1; % blur the posterior map
planeProb = bsxfun(@rdivide, planeProb, sum(planeProb, 3));

modelPlane.postProb = planeProb;

end

function vLine = vLineFromTwoVP(vp1, vp2)

A = cat(1, vp1, vp2);

[U S V] = svd(A, 0);
vLine = V(:,end);
vLine = vLine/vLine(3); % [h7, h8, 1]

end

repo: SelfExSR


text:
function test_bug2990

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

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug2990/4roboos_part1.mat'));
% the mri is from the dicom files, mri_realigned is in CTF coordinates, mri_realigned_resliced is after reslicing

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug2990/4roboos_part2.mat'));
% the cfg contains cfg.mri, which is more or less the mri_realigned (but more accurately coregistered than the sloppy ones in part1)
% the cfg also contains cfg.grid

cfg.checkconfig = 'loose';

cfg1 = cfg;
cfg1.mri = mri_realigned;
grid1 = ft_prepare_sourcemodel(cfg1);

cfg2 = cfg;
cfg2.mri = mri_realigned_resliced;
grid2 = ft_prepare_sourcemodel(cfg2);

figure; ft_plot_mesh(cfg.grid.template.pos(cfg.grid.template.inside,:))
figure; ft_plot_mesh(grid1.pos(grid1.inside,:))
figure; ft_plot_mesh(grid2.pos(grid2.inside,:))

% cfg.sourcemodel.template.pos(1,:)
% ans =
%     -8   -11    -7
%
% grid1.pos(1,:)
% ans =
%    -5.3781    8.3252   -0.2330
%
% Hmm, the first is in MNI, the second in CTF coordinates. Both are left-anterior-inferior.

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function pathway=markPathwayWithExpression(pathway,model,experiment,experimentOrder)
% markPathwayWithExpression
%   Adds gene expression data to a pathway object
%
%   pathway           map structure as generated by constructPathwayFromCelldesigner
%   model             a model structure
%   experiment        experiment structure as generated by getExpressionStructure
%   experimentOrder   cell array with two experiment names as defined in
%                     'experiment'
%
%   pathway           updated pathway object
%
%   Usage: pathway=markPathwayWithExpression(pathway,model,experiment,experimentOrder)

if numel(experimentOrder)~=2
    EM='This can only be done for two cases at the moment (experimentOrder must be two elements)';
    dispEM(EM);
end

%Check that experiment fit with experimentOrder
[present, expIds]=ismember(experimentOrder,experiment.experiments);

if ~all(present)
    EM='Not all experiments could be found in the experiment structure';
    dispEM(EM);
end

experiment.data=experiment.data(:,expIds);
experiment.experiments=experiment.experiments(:,expIds);

%Go through each species in pathway and if it's a reaction, get it's genes
%and log2-fold change in expression from model and experiment
for i=1:numel(pathway.listOfSpecies)
    if strcmp('PROTEIN', pathway.listOfSpecies(i).type)
        if isfield(pathway.listOfSpecies(i),'note')
            if ~isempty(pathway.listOfSpecies(i).note)
                %Get the reaction if present in model
                [present, index]=ismember(pathway.listOfSpecies(i).note,model.rxns);

                %If present, then get the genes
                if any(present)
                    [~, genes]=find(model.rxnGeneMat(index,:));

                    %If it was associated with genes match them to the ORFs
                    if any(genes)
                        [present, experimentIndexes]=ismember(model.genes(genes),experiment.orfs);

                        %Add annotation to pathway structure
                        if any(present)
                            pathway.listOfSpecies(i).orfs=experiment.orfs(experimentIndexes(present));
                            pathway.listOfSpecies(i).expA=experiment.data(experimentIndexes(present),1);
                            pathway.listOfSpecies(i).expB=experiment.data(experimentIndexes(present),2);
                        end
                    end
                end
            end
        end
    end
end
end

author: SysBioChalmers

repo: RAVEN


text:
function c = c8mat_add ( m, n, alpha, a, beta, b )

%*****************************************************************************80
%
%% C8MAT_ADD combines two C8MAT's with scalar factors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the order of A.
%
%    Input, complex ALPHA, the first scale factor.
%
%    Input, complex A(M,N), the first matrix.
%
%    Input, complex BETA, the second scale factor.
%
%    Input, complex B(M,N), the second matrix.
%
%    Output, complex C(M,N), the result.
%
  c(1:m,1:n) = alpha * a(1:m,1:n) + beta * b(1:m,1:n);

  return
end


text:
function line_ncc_rule_test ( )

%*****************************************************************************80
%
%% LINE_NCC_RULE_TEST tests the LINE_NCC_RULE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_NCC_RULE_TEST\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '  Test the LINE_NCC_RULE library.\n' );

  line_ncc_rule_test01 ( );
  line_ncc_rule_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_NCC_RULE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function [] = system_name1(block)

% Callback processing function for mask parameter
% - number of <Requirement name> used in the System
set_param(block, 'LinkStatus', 'none');
load_system('Simulink'); % Load Simulink library (background)
load_system('MIL_Test');

No = str2double(get_param(block,'No')); % Get # of VF-s
OldNo = length(find_system(block, 'regexp', 'on','LookUnderMasks','all','Name','<Requirement name>'));

for i = 1:max(No,OldNo)
    if OldNo > No && i > No && i <= OldNo % Remove
        % Remove lines
        set_param([block '/Bus Creator'], 'Inputs', num2str(No));
        SigAtTLines = get_param([block '/<Requirement name>' num2str(i)],'LineHandles');
        delete_line(SigAtTLines.Inport(1));
        delete_line(SigAtTLines.Outport(1));
        delete_line([block '/Arbitration'],['In' num2str(i) '/' num2str(1)],['Collect Verdicts/' num2str(i)]);
        % Remove blocks
        delete_block([block '/<Requirement name>' num2str(i)]);

        % remove lines from Bus Selector
        SigFrBS = get_param([block '/Bus Selector'], 'LineHandles');
        delete_line(SigFrBS.Outport(i));

        % remove Outport of Arbitration
        delete_block([block '/Arbitration/In' num2str(i)]);

        % Remove Inport of Bus Creator4

    elseif No > OldNo && i > OldNo && i <= No % Add
        % Add blocks
        add_block('MIL_Test/Validation Function/Test Evaluation Architecture/<system name>/<Requirement name>',[block '/<Requirement name>' num2str(i)]);
        set_param([block '/<Requirement name>' num2str(i)], 'LinkStatus', 'none');
        if i == 2
             pos = get_param([block '/<Requirement name>'], 'Position');
        else pos = get_param([block '/<Requirement name>' num2str(i-1)], 'Position');
        end
        pos1 = pos(1);
        pos3 = pos(3);
        pos2 = pos(2)+80;
        pos4 = pos(4)+80;
        set_param([block '/<Requirement name>' num2str(i)], 'Position',[pos1 pos2 pos3 pos4]);
        set_param([block '/Bus Creator'], 'Inputs', num2str(i));

        % Add Inports of Bus Creator4
        set_param([block '/Bus Creator4'], 'Inputs', num2str(i));

        % Add Outports of Bus Selector
        % OutSignal of Bus Seletor
        %Signal1 = 'signal1.signal1.signal1';
        %for j = 2 : i
        SignalP = strcat(',','signal',num2str(i),'.signal1.signal1');
        %Signal1 = strcat(Signal1,SignalP);
        Signal1 = get_param([block '/Bus Selector'], 'OutputSignals');
        Signal1 = strcat(Signal1,SignalP);
        %end
        set_param([block '/Bus Selector'], 'OutputSignals',Signal1)

        % Add lines
        add_line(block,['<Requirement name>' num2str(i) '/1'],['Bus Creator/' num2str(i)],'autorouting','on');
        add_line(block,['Signal Conversion/1'],['<Requirement name>' num2str(i) '/1'],'autorouting','on');

        % Add lines between Bus Selector and Bus Creator4
        add_line(block,['Bus Selector/' num2str(i)],['Bus Creator4/' num2str(i)],'autorouting','on');

        % Add Input to Arbitration
        add_block('Simulink/Ports & Subsystems/In1',[block '/Arbitration/In' num2str(i)]);
        if i == 2
             posi = get_param([block '/Arbitration/In1'], 'Position');
        else posi = get_param([block '/Arbitration/In' num2str(i-1)], 'Position');
        end
        pos1 = posi(1);
        pos3 = posi(3);
        pos2 = posi(2)+40;
        pos4 = posi(4)+40;
        set_param([block '/Arbitration/In' num2str(i)], 'Position',[pos1 pos2 pos3 pos4]);

        % fix position of Memory
        fix_pos([block '/Arbitration/Collect Verdicts'],[block '/Arbitration/Memory'],0,10);

        % Add Line from Bus Selector to Arbitration
        add_line(block,['Bus Selector/' num2str(i)],['Arbitration/' num2str(i)],'autorouting','on');

        % Set Inputports of 'Collect Verdicts'
        set_param([block '/Arbitration/Collect Verdicts'], 'Inputs', num2str(i+1));

        % Put the Memory Outsignal always in the last Inport of Collect
        % Verdicts
        delete_line([block '/Arbitration'],['Memory/' num2str(1)],['Collect Verdicts/' num2str(i)]);
        add_line([block '/Arbitration'],['Memory/' num2str(1)],['Collect Verdicts/' num2str(i+1)],'autorouting','on');

        % Add line from 'in' to 'Collect Verdicts', which happens in the
        % undersystem Arbitration
        add_line([block '/Arbitration'],['In' num2str(i) '/' num2str(1)],['Collect Verdicts/' num2str(i)],'autorouting','on');
    end
end

% Remove Outport of Bus Selector & Input of Bus Creator4

    if OldNo > No  % Remove
        set_param([block '/Bus Creator4'], 'Inputs', num2str(No));
        delete_line([block '/Arbitration'],['Memory/' num2str(1)],['Collect Verdicts/' num2str(OldNo+1)]);
        set_param([block '/Arbitration/Collect Verdicts'], 'Inputs', num2str(No+1));
        add_line([block '/Arbitration'],['Memory/' num2str(1)],['Collect Verdicts/' num2str(No+1)],'autorouting','on');
        if No == 1
            S1 = 'signal1.signal1.signal1';
        else
            S1 = 'signal1.signal1.signal1';
            for j = 2:No
              S1 = strcat(S1,',','signal',num2str(j),'.signal1.signal1');
           end
        end
        set_param([block '/Bus Selector'], 'OutputSignals',S1);

    end


text:
function [q]=zn2fpd(input)
% [q]=zn2fpd(input)
% Ziegler-Nichols controller for 2nd order processes
% This function computes parameters of the controller (q0, q1, q2, q3, q4).
% Controller is based on forward rectangular method of discretization
% replacing derivation by a four-point difference.
% Transfer function of the controller is as follows:
%
%            q0 + q1*z^-1 + q2*z^-2 + q3*z^-3 + q4*z^-4
% G(z^-1) = --------------------------------------------
%                           1 - z^-1
%
% Transfer function of the controlled system is:
%
%               b1*z^-1 + b2*z^-2
% Gs(z^-1) = -----------------------
%             1 + a1*z^-1 + a2*z^-2
%
% Input: input ... input parameters
%                  input(1) ... a1
%                  input(2) ... b1
%                  input(3) ... a2
%                  input(4) ... b2
%                  input(5) ... sample time T0
% Output: qp ... controller parameters
%                qp(1) ... q0
%                qp(2) ... q1
%                qp(3) ... q2
%                qp(4) ... p1
%                qp(5) ... p2

a1 = input(1);
b1 = input(2);
a2 = input(3);
b2 = input(4);
T0 = input(5);

% compute ultimate gain and frequency
[Kpu, Tu] =  ultim([b1 b2],[a1 a2],T0);

Kp = 0.6*Kpu;
Ti = Tu/2;
Td = Tu/8;

q0 = Kp*(1 + T0/Ti + Td/(6*T0));
q1 = Kp*(-1 + Td/(3*T0));
q2 = Kp*(-Td/T0);
q3 = Kp*(Td/(3*T0));
q4 = Kp*(Td/(6*T0));

q=[q0;q1;q2;q3;q4];


text:
function [sys,x0,str,ts]=Plant1_x1x2wt(t,x,u,flag)

switch flag,
    case 0,
        [sys,x0,str,ts]=mdlInitializeSizes;
    case 1,
        sys=mdlDerivatives(t,x,u);
    case 3,
        sys=mdlOutputs(x);
    case {2,4,9},
        sys=[];
    otherwise
        error(['Unhandled flag=',num2str(flag)]);
end
function [sys,x0,str,ts]=mdlInitializeSizes
    sizes=simsizes;
    sizes.NumContStates=2;
    sizes.NumDiscStates=0;
    sizes.NumOutputs=1;
    sizes.NumInputs=1;
    sizes.DirFeedthrough=1;
    sizes.NumSampleTimes=1;
    sys=simsizes(sizes);
    x0=[0;0];
    str=[];
    ts=[0 0];
function sys=mdlDerivatives(t,x,u)
    sys(1)=x(2);
    sys(2)=-5*x(1)-20*x(2)++0.5*sign(sin(t))+u;
function sys=mdlOutputs(x)
    sys=x(1);

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
% op_alignrcvrs.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% [out,coilcombos]=op_alignrcvrs(in,point,mode,coilcombos);
%
% DESCRIPTION:
% phase align the receiver channels without combining them.
%
% INPTUS:
% in            = input spectrum in matlab structure format.
% point         = Index of point in time domain to use for phase reference.
%                 (optional.  Default = 1);
% mode          = Method for estimating the coil weights and phases (optional.  Default = 'w').
%                 -'w' performs amplitude weighting of channels based on the
%                 maximum signal of each coil channel.
%                 -'h' performs amplitude weighting of channels based on the
%                 maximum signal of each coil channel divided by the square of
%                 the noise in each coil channel (as described by Hall et al.
%                 Neuroimage 2014).
% coilcombos	= (optional)  The predetermined coil phases and amplitudes as
%                 generated by the op_getcoilcombos.m function.  If this
%                 argument is provided, the 'point', and 'mode', arguments
%                 will be ignored.
%
% OUTPUTS:
% out           = Output following alignment of rf channels.
% coilcombos    = Structure containing two fields:
%                   ph:  Vector of coil phases (in degrees) used for alignment.
%                   sig: Vector of coil weights.

function [out,coilcombos]=op_alignrcvrs(in,point,mode,coilcombos);

if in.flags.addedrcvrs
    error('ERROR:  Receivers have already been combined!  Aborting!');
end

%To get best possible SNR, add the averages together (if it hasn't already been done):
if in.dims.averages>0
    av=op_averaging(in);
else
    av=in;
end

%also, for best results, we will combine all subspectra:
if nargin<4
    if in.flags.isFourSteps
        av=op_fourStepCombine(av);
    end
    if in.dims.subSpecs>0
        av=op_combinesubspecs(av,'summ');
    end
    if nargin < 3
        mode='w';
        if nargin < 2
            point=1;
        end
    end

end
avfids=av.fids;
avspecs=av.specs;

%initialize phase matrix and the amplitude maxtrix that are the size of nPoints x Coils
ph=ones(in.sz(in.dims.t),in.sz(in.dims.coils));
sig=ones(in.sz(in.dims.t),in.sz(in.dims.coils));

if nargin<4
    coilcombos.ph=zeros(in.sz(in.dims.coils),1);
    coilcombos.sig=zeros(in.sz(in.dims.coils),1);
end

%now start finding the relative phases between the channels and populate
%the ph matrix
for n=1:in.sz(in.dims.coils)
    if nargin<4
        p=phase(avfids(point,n,1,1));
        coilcombos.ph(n)=p;
        ph(:,n)=p*ph(:,n);
        switch mode
            case 'w'
                S=abs(avfids(point,n,1,1));
                coilcombos.sig(n)=S;
                sig(:,n)=S*sig(:,n);
            case 'h'
                S=abs(avfids(point,n,1,1));
                N=std(avfids(end-100:end,n,1,1));
                sig(:,n)=(S/(N.^2))*sig(:,n);
                coilcombos.sig(n)=(S/(N.^2));
        end
    else
        ph(:,n)=coilcombos.ph(n)*ph(:,n);
        sig(:,n)=coilcombos.sig(n)*sig(:,n);
    end
end

%now replicate the phase matrix to equal the size of the original matrix:
replicate=in.sz;
replicate(1)=1;
replicate(2)=1;
ph=repmat(ph,replicate);
%sig=repmat(sig,replicate);
%sig=sig/max(max(max(max(sig))));

%now apply the phases by multiplying the data by exp(-i*ph);
fids=in.fids.*exp(-i*ph);
fids_presum=fids;
specs_presum=fftshift(ifft(fids,[],in.dims.t),in.dims.t);

%FILLING IN DATA STRUCTURE
out=in;
out.fids=fids_presum;
out.specs=specs_presum;

author: CIC-methods

repo: FID-A

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


text:
function [l, L_rf, L_sf, L_k, L_seg, L_n] = ...
    retroProjIdpLinFromPinHoleOnRob(Rf, Sf, k, seg, n)

% RETROPROJIDPLINFROMPINHOLEONROB retroprj Idp Line from pinhole on robot.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if nargout == 1

    ls = invPinHoleIdpLin(k, seg, n) ;
    lr = fromFrameIdpLin(Sf, ls);
    l  = fromFrameIdpLin(Rf, lr);

else % Jacobians requested

    [ls, LS_seg, LS_n, LS_k] = invPinHoleIdpLin(seg, n, k) ;
    [lr, LR_sf, LR_ls]       = fromFrameIdpLin(Sf, ls);
    [l, L_rf, L_lr]          = fromFrameIdpLin(Rf, lr);

    L_sf  = L_lr*LR_sf;
    L_ls  = L_lr*LR_ls;
    L_k   = L_ls*LS_k;
    L_seg = L_ls*LS_seg;
    L_n   = L_ls*LS_n;

end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function ho=addmarkers(hi,n)
% ADDMARKERS    Add a pre-defined number of markers to a plot instead of
% one at each datapoint which is default in Matlab.
%
% Usage:
% ADDMARKERS(HI) adds 5 markers to the input handle HI, J=1..length(HI)
% ADDMARKERS(HI,N) adds N markers to the input handle HI, J=1..length(HI)
% HO=ADDMARKERS(...) returns the output handle HO for future usage such as
% changing colors, markers, etc.
%
% Hint: use the output handle when adding legends, e.g.
% LEGEND(HO,'legend 1',...,'legend N')
%
% The markers will be added in the following order:
%
% o     Circle
% s     Square
% d     Diamond
% x     Cross
% +     Plus sign
% *     Asterisk
% .     Point
% ^     Upward-pointing triangle
% v     Downward-pointing triangle
% >     Right-pointing triangle
% <     Left-pointing triangle
% p     Five-pointed star (pentagram)
% h     Six-pointed start (hexagram)
%
% By Matthijs Klomp at Saab Automobile AB, 2009-10-11
% E-mail: matthijs.klomp@gm.com
%
if nargin<2, n=5; end % default number of markers
if nargin<1, error('Supply an input handle as input argument.'), end
figure(get(get(hi(1),'parent'),'parent'))% plot in the figure of the handle
subplot(get(hi(1),'parent'))            % plot in the subplot of the handle
hold on                                 % do not overwrite the current plot
markers = {'o','s','d','^','v','x','+','*','.','>','<','p','h'};
ho = hi;                                % initialize output handle
for i = 1:length(hi)
    x = get(hi(i),'xdata');             % get the independent variable
    y = get(hi(i),'ydata');             % get the dependent variable data
    s = linspace(1,n,length(x));        % sampling independent variable
    sn = [1 (2:n-1)+randn(1,n-2)/n n];  % add some noise to avoid overlap
    xrs = interp1(s,x,sn,'nearest');    % downsample to n datapoints
    yrs = interp1(s,y,sn,'nearest');    % downsample to n datapoints
    plot(xrs,yrs,...                    % Plot the markers
        'Marker',markers{i},'LineStyle','None','Color',get(hi(i),'Color'));
    ho(i) = plot([0 1],[1 1]*NaN,...    % Create the output handle
        'Marker',markers{i},...
        'LineStyle',get(hi(i),'Linestyle'),...
        'Color',get(hi(i),'Color'));
end
if nargout==0, clear ho, end


text:
function [DataMat, ChannelMat] = in_data_ws_csv(DataFile)
% IN_DATA_WS_CSV: Imports a Wearable Sensing CSV file.

% @=============================================================================
% 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, 2019

% ===== READ CSV =====
bst_progress('text', 'Reading csv file...');
% Open file
fid = fopen(DataFile, 'r');
if (fid == -1)
    error('Cannot open file.');
end
% Read header line with channel names
Labels = [];
while 1
    csvLine = fgetl(fid);
    % End of file
    if (csvLine == -1)
        break;
    % Empty line
    elseif isempty(csvLine)
        continue;
    % Comment line
    elseif (csvLine(1) == '#')
        continue;
    % Header line with channel names
    elseif (nnz(csvLine == ',') > 0)
        Labels = strtrim(str_split(csvLine, ','));
        break;
    end
end
% Channel names not found
if isempty(Labels)
    fclose(fid);
    error('Header line with channel names not found in CSV.');
end
% Default reading pattern
cellFormat = repmat({'%f'}, 1, length(Labels));

% Final reading format for the rest of the file
strFormat = sprintf('%s ', cellFormat{:});
strFormat(end) = [];
% Read the rest of the file
RecMat = textscan(fid, strFormat, 'Delimiter', ',');
% Close file
fclose(fid);
% Check for errors
if isempty(RecMat) || isempty(RecMat{1})
    error('File could not be read as CSV.');
end
% If the last line is incomplete: delete it
lastFullRow = length(RecMat{end});
if (length(RecMat{1}) > lastFullRow)
    for i = 1:length(RecMat)
        if (length(RecMat{i}) > lastFullRow)
            RecMat{i} = RecMat{i}(1:lastFullRow);
        end
    end
end

% ===== GET TIME =====
bst_progress('text', 'Processing time stamps...');
% Find timestamp
iTimestamp = find(strcmpi(Labels, 'Time'));
if isempty(iTimestamp)
    error('Missing time column.');
end
% Convert time stamps to datenum
rawTime = RecMat{iTimestamp}';
% Set t=0 at the first sample
rawTime = rawTime - rawTime(1);
% Detect sampling frequency
sfreq = 1./mean(diff(rawTime));
sfreq = round(sfreq * 100) ./ 100;
% Reconstructed file time
Time = (0:length(rawTime)-1) ./ sfreq;

% ===== GET DATA =====
bst_progress('text', 'Reconstructing data matrix...');
% Get other useless channels
iComments = find(strcmpi(Labels, 'Comments'));
% Find a column that has non-empty values that is not a special column
iColAll = setdiff(1:length(Labels), [iTimestamp, iComments]);
iColRec = iColAll(~cellfun(@(c)all(isnan(c)), RecMat(iColAll)));
% Rebuild data matrix
rawF = [RecMat{iColRec}]';
% Remove other missing values: replacing with values
iMissing = find(any(isnan(rawF),1));
if ~isempty(iMissing)
    disp(sprintf('BST> Muse: Missing data at %d time points. Replacing with zeros...', length(iMissing)));
    rawF(isnan(rawF)) = 0;
end
nChannels = size(rawF,1);

% % ===== DETECT EVENTS =====
% % Find events
% iTimeEvt = find(all(isnan(rawF),1));
% % Get event timing
% evtTime = rawTime(iTimeEvt);
% % Get event labels
% if ~isempty(iElements) && ~all(cellfun(@isempty, RecMat{iElements}(iTimeEvt)))
%     evtLabels = cellfun(@(c)strrep(c, '/muse/elements/', ''), RecMat{iElements}(iTimeEvt), 'UniformOutput', 0);
% else
%     evtLabels = repmat({'Trigger'}, length(iTimeEvt), 1);
% end
% Remove events from data matrix
% rawF(:,iTimeEvt) = [];
% rawTime(iTimeEvt) = [];

% ===== BRAINSTORM DATA =====
% Get file name
[fPath, fBase, fExt] = bst_fileparts(DataFile);
% Create empty structure
DataMat = db_template('DataMat');
% Fill structure
DataMat.F           = rawF .* 1e-6;  % Data stored in microV
DataMat.Time        = Time;
DataMat.Comment     = fBase;
DataMat.ChannelFlag = ones(nChannels, 1);
DataMat.nAvg        = 1;
DataMat.Device      = 'Wearable Sensing';
% % Add events
% if ~isempty(iTimeEvt)
%     % Initialize events list
%     DataMat.Events = repmat(db_template('event'), 0);
%     % Events list
%     uniqueEvt = unique(evtLabels);
%     % Build events list
%     for iEvt = 1:length(uniqueEvt)
%         % Find all the occurrences of this event
%         iOcc = find(strcmpi(uniqueEvt{iEvt}, evtLabels));
%         % Set event
%         DataMat.Events(iEvt).label    = strtrim(uniqueEvt{iEvt});
%         DataMat.Events(iEvt).times    = unique(round(evtTime(iOcc) .* sfreq)) ./ sfreq;
%         DataMat.Events(iEvt).epochs   = 1 + 0*DataMat.Events(iEvt).times;
%         DataMat.Events(iEvt).select   = 1;
%         DataMat.Events(iEvt).channels = [];
%         DataMat.Events(iEvt).notes    = [];
%     end
% end

% ===== CHANNEL FILE =====
ChannelMat = db_template('channelmat');
for i = 1:nChannels
    % Split type/name
    ch = Labels{iColAll(i)};
    ChannelMat.Channel(i).Name = ch;
    switch (ch)
        case 'Trigger'
            ChannelMat.Channel(i).Type = 'Trigger';
        case 'Event'
            ChannelMat.Channel(i).Type = 'Event';
        case 'Time_Offset'
            ChannelMat.Channel(i).Type = 'Time_Offset';
        case 'CM'
            ChannelMat.Channel(i).Type = 'CM';
        case {'ADC_Status', 'ADC_Sequence'}
            ChannelMat.Channel(i).Type = 'Trigger';
        otherwise
            ChannelMat.Channel(i).Type = 'EEG';
    end
    ChannelMat.Channel(i).Loc     = [];
    ChannelMat.Channel(i).Orient  = [];
    ChannelMat.Channel(i).Comment = '';
    ChannelMat.Channel(i).Weight  = 1;
end

author: brainstorm-tools


text:
function [data,eid,name] = DBGetValues(query,direction)

%DBGetValues - Group values for all variables that match given criteria.
%
%  Concatenate individual values for all variables matching the optional query.
%  Values must be scalars, vectors or matrices. Concatenation can be performed
%  horizontally or vertically. Matrices can also be concatenated along the
%  third dimension (z). By default, values of size MxN are concatenated
%  vertically if M<=N, or horizontally otherwise.
%
%  USAGE
%
%    [data,eid,name] = DBGetValues(query,direction)
%
%    query          optional list query (WHERE clause; see Example)
%    direction      optional concatenation direction ('h', 'v' or 'z')
%
%  OUTPUT
%
%    data           concatenation of individual values
%    eid            experiment ID (identifier string)
%    name           descriptive name (identifier string)
%
%  EXAMPLE
%
%    In this example, place cells were recorded on successive days while the
%    animal explored a maze, then during sleep. Firing fields were stored in a
%    database using eids like '20120213-Maze-(1,2)' (where 1,2 corresponds to
%    tetrode 1, cluster 2) and named 'FiringFields'. Mean firing rates during
%    sleep were named 'MeanRate'.
%
%    Get all firing fields (for details on query syntax, see an SQL manual):
%
%    [fields,eid1] = DBGetValues('eid like "%Maze%" and name="FiringField"');
%
%    Get all firing rates during sleep:
%
%    [rates,eid2] = DBGetValues('eid like "%Sleep%" and name="MeanRate"');
%
%  SEE
%
%    See also DBMatchValues, DBGetVariables, DBAddVariable, DBGetFigures.
%

% Copyright (C) 2007-2013 by Michaël Zugaro
%
% 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.

% Make sure MyM is installed and functional
CheckMyM;

data = [];
eid = {};
name = {};

% Parameters
if nargin < 1,
	query = '';
	direction = [];
elseif nargin == 1,
	if strcmp(lower(query),'h'),
		query = '';
		direction = 'h';
	elseif strcmp(lower(query),'v'),
		query = '';
		direction = 'v';
	elseif strcmp(lower(query),'z'),
		query = '';
		direction = 'z';
	else
		direction = [];
	end
end

% Query database
results = DBGetVariables(query,'output','full');
if isempty(results.v), return; end
eid = results.eid;
name = results.name;

first = results.v{1};
if ~isnumeric(first),
	error('Values are not numerical and cannot be concatenated (type ''help DBGetValues'' for details).');
end

% Automatic concatenation direction
if isempty(direction),
	if size(first,1) <= size(first,2),
		direction = 'v';
	else
		direction = 'h';
	end
end

% Concatenation

% Determine sizes of retrieved values (arrays)
sizes = cellfun(@size,results.v,'UniformOutput',false);
% Make sure they all have the same number of dimensions
nDims = cellfun(@length,sizes);
if any(nDims~=nDims(1)),
	error('All values should have the same number of dimensions (type ''help DBGetValues'' for details).');
end
% Determine 'circumscribed' array size
s = max(unique(vertcat(sizes{:}),'rows'),[],1);
% Extend all arrays
results.v = cellfun(@(x) ExtendArray(x,s),results.v,'uniformoutput',false);
% Concatenate in the appropriate direction
if strcmp(direction,'v'),
	data = vertcat(results.v{:});
elseif strcmp(direction,'h'),
	data = [results.v{:}];
else
	data = cat(3,results.v{:});
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function allw = traceToWaveform(obj)
   %traceToWaveform converts traceobjects to waveforms
   %    w = traceToWaveform(antelopesource)
   %
   % Note: this may return multiple waveform objects, depending upon
   % how many segments and/or scnl's.

   wBlank = antelopesource.blankWaveformWithCalib();

   try
      traceCount = dbnrecs(obj.trpointer);
      assert(traceCount > 0);
   catch er
      disp('traceToWaveform: No traces found to convert');
      allw = wBlank([]);
      return
   end

   % preallocations
   spikes(traceCount).mask = false; %used to track data spikes & inf values
   allw = repmat(wBlank,traceCount,1);

   % LOOP twice through segments represented by this trace object
   % 1st loop: find signal spikes, and get header info.
   for seg = 1:traceCount
      obj.trpointer.record = seg - 1;
      allw(seg) = waveformMetadata(wBlank, obj.trpointer);
      spikes(seg).mask = findSpikes(obj.trpointer);
   end

   allw = applyCalibrationsToAll(trpointer, allw);

   % 2nd loop: assign data to the waveforms.
   replaceWithNan = @(W,BAD) setsamples(W, BAD.mask, nan);
   for seg = 1:traceCount
      obj.trpointer.record = seg - 1;
      allw(seg) = set(allw(seg), 'data', trextract_data(obj.trpointer));
      allw(seg) = replaceWithNan(allw(seg), spikes(seg));
   end
end

function w = waveformMetadata(w, trpointer)
   s = db2struct(trpointer); %do once, get one for each segment
   [sunit, ~] = segtype2units(s.segtype);
   w = set(w, ...
      'station', s.sta, ...
      'channel', s.chan, ...
      'network', s.net, ...
      ...'location', s.loc, ... % unknown if 'loc' really is a field
      'start', datenum(epoch2str(s.time,'%Y-%m-%d %H:%M:%S.%s')), ...
      'freq', s.samprate, ...
      'units', sunit, ...
      'calib', s.calib);
end

function spikes = findSpikes(trpointer)
      % data spikes must be known PRIOR to applying calibration
      maxAllowableSignal = (realmax('single') * 1e-2);
      data = trextract_data(trpointer);
      spikes = (abs(data) >= maxAllowableSignal) | isinf(data);
end

function allw = applyCalibrationsToAll(trpointer, allw)
   % now, apply calibrations to all traces at once
   trapply_calib(trpointer);
   hasCalibs = get(allw,'calib') ~= 0;
   allw(hasCalibs) = set(allw(hasCalibs),'calibration_applied','YES');
end

author: geoscience-community-codes

repo: GISMO


text:
%% Find Image Rotation and Scale Using Automated Feature Matching
% This example shows how to automatically align two images that differ by a
% rotation and a scale change. It closely parallels another example titled
% <matlab:showdemo('RotationFitgeotransExample') Find Image Rotation and Scale>.
% Instead of using a manual approach to register the two images, it
% utilizes feature-based techniques found in the Computer Vision System
% Toolbox(TM) to automate the registration process.
%
% In this example, you will use |detectSURFFeatures| and
% |vision.GeometricTransformEstimator| System object to recover rotation
% angle and scale factor of a distorted image. You will then transform the
% distorted image to recover the original image.

% Copyright 1993-2014 The MathWorks, Inc.

function recovered = homographyAlignment(original, distorted, verbose)
% clc; clear; close all

%% Step 1: Read Image
% Bring an image into the workspace.
if isempty(original)
    original = imread('00001.jpg');
end
original_c = original;
original = rgb2gray(original);
if verbose
    imshow(original);
    text(size(original,2),size(original,1)+15, ...
        'Image courtesy of Massachusetts Institute of Technology', ...
        'FontSize',7,'HorizontalAlignment','right');
end

%% Step 2: Resize and Rotate the Image
if isempty(distorted)
    % scale = 0.7;
    % J = imresize(original, scale); % Try varying the scale factor.
    %
    % theta = 30;
    % distorted = imrotate(J,theta); % Try varying the angle, theta.
    distorted = rgb2gray(imread('00003.jpg'));
end
distorted_c = distorted;
distorted = rgb2gray(distorted);
if verbose
    figure, imshow(distorted)
end

%%
% You can experiment by varying the scale and rotation of the input image.
% However, note that there is a limit to the amount you can vary the scale
% before the feature detector fails to find enough features.

%% Step 3: Find Matching Features Between Images
% Detect features in both images.
ptsOriginal  = detectSURFFeatures(original);
ptsDistorted = detectSURFFeatures(distorted);

%%
% Extract feature descriptors.
[featuresOriginal,   validPtsOriginal]  = extractFeatures(original,  ptsOriginal);
[featuresDistorted, validPtsDistorted]  = extractFeatures(distorted, ptsDistorted);

%%
% Match features by using their descriptors.
indexPairs = matchFeatures(featuresOriginal, featuresDistorted);

%%
% Retrieve locations of corresponding points for each image.
matchedOriginal  = validPtsOriginal(indexPairs(:,1));
matchedDistorted = validPtsDistorted(indexPairs(:,2));

%%
% Show point matches. Notice the presence of outliers.
if verbose
    figure;
    showMatchedFeatures(original,distorted,matchedOriginal,matchedDistorted);
    title('Putatively matched points (including outliers)');
end

%% Step 4: Estimate Transformation
% Find a transformation corresponding to the matching point pairs using the
% statistically robust M-estimator SAmple Consensus (MSAC) algorithm, which
% is a variant of the RANSAC algorithm. It removes outliers while computing
% the transformation matrix. You may see varying results of the
% transformation computation because of the random sampling employed by the
% MSAC algorithm.
[tform, inlierDistorted, inlierOriginal, status] = estimateGeometricTransform(...
    matchedDistorted, matchedOriginal, 'projective');
if status
    recovered = original_c;
    return;
end

%%
% Display matching point pairs used in the computation of the
% transformation matrix.
if verbose
    figure;
    showMatchedFeatures(original,distorted, inlierOriginal, inlierDistorted);
    title('Matching points (inliers only)');
    legend('ptsOriginal','ptsDistorted');
end

%% Step 5: Solve for Scale and Angle
% Use the geometric transform, TFORM, to recover
% the scale and angle. Since we computed the transformation from the
% distorted to the original image, we need to compute its inverse to
% recover the distortion.
%
%  Let sc = scale*cos(theta)
%  Let ss = scale*sin(theta)
%
%  Then, Tinv = [sc -ss  0;
%                ss  sc  0;
%                tx  ty  1]
%
%  where tx and ty are x and y translations, respectively.
%

%%
% Compute the inverse transformation matrix.
Tinv  = tform.invert.T;

ss = Tinv(2,1);
sc = Tinv(1,1);
scale_recovered = sqrt(ss*ss + sc*sc);
theta_recovered = atan2(ss,sc)*180/pi;

%%
% The recovered values should match your scale and angle values selected in
% *Step 2: Resize and Rotate the Image*.

%% Step 6: Recover the Original Image
% Recover the original image by transforming the distorted image.
outputView = imref2d(size(original_c));
recovered = zeros(size(original_c));

[h,w,l] = size(original_c);
[mx,my] = meshgrid(1:w,1:h);
tmp = [mx(:)';my(:)';ones(1,h*w)];
tmp = (Tinv)'*tmp;
X = reshape(tmp(1,:)./tmp(3,:),h,w);
Y = reshape(tmp(2,:)./tmp(3,:),h,w);
recovered = ba_interp2(single(distorted_c),X,Y,'linear');

% for i=1:l
%     recovered(:,:,i) = interp2(mx,my,distorted_c(:,:,i),X,Y,'linear',0);
% end

% for ic = 1:size(original_c,3)
%     recovered(:,:,ic)  = imwarp(squeeze(distorted_c(:,:,ic)),tform,'OutputView',outputView,'FillValues',0,'SmoothEdges',true);
% end

%%
% Compare |recovered| to |original| by looking at them side-by-side in a montage.
if verbose
    figure, imshowpair(original_c,recovered,'montage')
%%
% The |recovered| (right) image quality does not match the |original| (left)
% image because of the distortion and recovery process. In particular, the
% image shrinking causes loss of information. The artifacts around the edges are
% due to the limited accuracy of the transformation. If you were to detect
% more points in *Step 4: Find Matching Features Between Images*,
% the transformation would be more accurate. For example, we could have
% used a corner detector, |detectFASTFeatures|, to complement the SURF
% feature detector which finds blobs. Image content and image size also
% impact the number of detected features.

    displayEndOfDemoMessage(mfilename)
end

author: shuochsu

repo: DeepVideoDeblurring


text:
function s = ymdf_to_s_islamic ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_S_ISLAMIC writes an Islamic YMDF date into a string.
%
%  Format:
%
%    AH YYYY/MM/DD.FF
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, string S, a representation of the date.
%
  s = sprintf ( 'AH %d/%02d/%05.2f', y, m, d + f );

  return
end


text:
clear; close all; clc;

n = 150;
xx = linspace(-4,4,n);
my_polynomial = [1, -2, -3, 4];
yy = polyval(my_polynomial, xx);

n_step = 150;

newVid = VideoWriter('pic1', 'MPEG-4'); % New
newVid.FrameRate = 30;
newVid.Quality = 100;
open(newVid);

figure('position',[490, 110, 900, 650],'color','w');
for i_step = 1:n_step
    yy2plot = (yy - zeros(1, n)) * i_step / n_step + zeros(1, n);
    plot2DPlane(-4,4,-6,6,false);
    hold on;

    if i_step == 1

        plot(xx, zeros(1,n),'linewidth',2);
        xlim([-3, 4])
        ylim([-6, 6])
        grid on;
        xlabel('$$x$$','interpreter','latex');
        ylabel('$$y = f(x)$$','interpreter','latex');
        title('$$y=f(x)=x^3-2x^2-3x+4$$','interpreter','latex')

        for i = 1:29 % 첫 장면에서 1.5초 대기할 수 있도록

            writeVideo(newVid, getframe(gcf))
        end
        drawnow;
        cla;
    end
    plot2DPlane(-4,4,-6,6,false);

    plot(xx, yy2plot,'linewidth',2);
    xlim([-3, 4])
    ylim([-6, 6])
    grid on;
    xlabel('$$x$$','interpreter','latex');
    ylabel('$$y = f(x)$$','interpreter','latex');
    title('$$y=f(x)=x^3-2x^2-3x+4$$','interpreter','latex')
    drawnow;

    writeVideo(newVid, getframe(gcf))

    if i_step < n_step
        cla;
    end
end

for i = 1:30 % 마지막 장면에서 1.5초 더 대기할 수 있도록
    writeVideo(newVid, getframe(gcf))
end

close(newVid)

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
% @author: Maziar Raissi

function KDV()
clc; close all;

plt = 1;
plt_pred = 0;
save_plt = 1;

addpath ..
addpath ../Utilities
addpath ../Kernels/KDV
addpath ../Utilities/export_fig

function CleanupFun()
    rmpath ..
    rmpath ../Utilities
    rmpath ../Kernels/KDV
    rmpath ../Utilities/export_fig
end

finishup = onCleanup(@() CleanupFun());

set(0,'defaulttextinterpreter','latex')

%% Load Data
load('../Data/kdv.mat', 'usol', 't', 'x')
u_star = real(usol); % 512x201
t_star = t; % 201x1
x_star = x';   % 512x1
N_star = size(x_star,1);
nsteps = size(t_star,1)-1;

%% Setup
N0 = 111;
N1 = 109;
%% Clean Data
rng('default')
i = randi(nsteps);
dt = t_star(i+1) - t_star(i);

idx0 = randsample(N_star, N0);
x0 = x_star(idx0,:);
u0 = u_star(idx0,i);

idx1 = randsample(N_star,N1);
x1 = x_star(idx1,:);
u1 = u_star(idx1,i+1);

hyp = [log([1.0 1.0]) 0.0 0.0 -4.0];
model = HPM(x1, u1, x0, u0, dt, hyp);
model = model.train(5000);

hyp = model.hyp;
params = hyp(3:4);

[pred_n_star, var_n_star] = model.predict(x_star);
var_n_star = abs(diag(var_n_star));

error = norm(pred_n_star - u_star(:,i+1))/norm(u_star(:,i+1));

fprintf(1,'=========================\n');
fprintf(1,'Step: %d, Time = %.2f\n\nNLML = %.2f, Error = %.2e\n\n', i, ...
    t_star(i+1), model.NLML, error);

str = sprintf('%.4f  ', params);
fprintf('Parameters: %s\n\n', str)
fprintf(1,'=========================\n\n');

if plt_pred == 1
    figure();
    plot_prediction_1D(x_star, u_star(:,i+1), pred_n_star, var_n_star, ...
        '$x$', '$u(t,x)$', 'Prediction (clean data)');

    drawnow;
end

%% Noisy Data
noise = 0.01;
u0 = u0 + noise*std(u0)*randn(size(u0));
u1 = u1 + noise*std(u1)*randn(size(u1));

hyp = [log([1.0 1.0]) 0.0 0.0 -4.0];
model = HPM(x1, u1, x0, u0, dt, hyp);
model = model.train(5000);

hyp = model.hyp;
params_noise = hyp(3:4);

[pred_n_star, var_n_star] = model.predict(x_star);
var_n_star = abs(diag(var_n_star));

error = norm(pred_n_star - u_star(:,i+1))/norm(u_star(:,i+1));

fprintf(1,'=========================\n');
fprintf(1,'Step: %d, Time = %.2f\n\nNLML = %.2f, Error = %.2e\n\n', i, ...
    t_star(i+1), model.NLML, error);

str = sprintf('%.4f  ', params_noise);
fprintf('Parameters: %s\n\n', str)
fprintf(1,'=========================\n\n');

if plt_pred == 1
    figure();
    plot_prediction_1D(x_star, u_star(:,i+1), pred_n_star, var_n_star, ...
        '$x$', '$u(t,x)$', 'Prediction (noisy data)');

    drawnow;
end

%% Plot Results

if plt == 1
    fig = figure();
    set(fig,'units','normalized','outerposition',[0 0 1 .5])
    subplot(3,2,1:2)
    plot_surface(t_star, x_star, u_star, '$t$', '$x$', '$u(t,x)$');

    hold on
    plot3([t_star(i) t_star(i)],get(gca,'ylim'),[10 10],'w','LineWidth',2)
    plot3([t_star(i+1) t_star(i+1)],get(gca,'ylim'),[10 10],'w','LineWidth',2)

    subplot(3,2,3);
    tit = sprintf('$t = $ %.2f\n%d training data\n', t_star(i), N0);
    plot_data_1D(x_star, u_star(:,i), x0, u0, '$x$', '$u(t,x)$', tit);

    subplot(3,2,4);
    tit = sprintf('$t = $ %.2f\n%d training data\n', t_star(i+1), N1);
    plot_data_1D(x_star, u_star(:,i+1), x1, u1, '$x$', '$u(t,x)$', tit);

    subplot(3,2,5:6);
    s1 = '$\begin{tabular}{ |c|c| }  \hline Correct PDE & $u_t + 6 u u_x + u_{xxx} = 0$ \\  \hline Identified PDE (clean data) &';
    s2 = sprintf('$u_t + %.3f u u_x + %.3f u_{xxx} = 0$', params(1), params(2));
    s3 = ' \\  \hline Identified PDE (1\% noise) &';
    s4 = sprintf('$u_t + %.3f u u_x + %.3f u_{xxx} = 0$', params_noise(1), params_noise(2));
    s5 = ' \\  \hline \end{tabular}$';
    s = strcat(s1,s2,s3,s4,s5);
    text(0.1,0.8,s,'interpreter','latex','FontSize',18)
    axis off

    if save_plt == 1
        export_fig ../Figures/KDV.png -r300
    end

    drawnow();
end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
function p = getExpParams(predidx)

switch predidx
    case 0
        p.name = 'Pishchulin et., CVPR''13';
        p.predFilename = '';
    case 1
        p.name = 'Pishchulin et., ICCV''13';
        p.predFilename = './pred/pishchulin13iccv/pred_keypoints_lsp_oc';
        p.colorIdxs = [1 1];
    case 2
        p.name = 'Tompson et al., NIPS''14';
        p.predFilename = './pred/tompson14nips/pred_keypoints_lsp_pc';
        p.colorIdxs = [2 1];
    case 3
        p.name = 'Chen&Yuille, NIPS''14';
        p.predFilename = './pred/chen14nips/pred_keypoints_lsp_oc';
        p.colorIdxs = [5 1];
    case 4
        p.name = 'Ramakrishna et al., ECCV''14';
        p.predFilename = './pred/ramakrishna14eccv/pred_keypoints_lsp_oc';
        p.colorIdxs = [4 1];
    case 5
        p.name = 'Ouyang et al., CVPR''14';
        p.predFilename = './pred/ouyang14cvpr/pred_keypoints_lsp_oc';
        p.colorIdxs = [7 1];
    case 6
        p.name = 'Pishchulin et., ICCV''13';
        p.predFilename = './pred/pishchulin13iccv/pred_keypoints_lsp_pc';
        p.colorIdxs = [1 1];
    case 7
        p.name = 'Pishchulin et., ICCV''13';
        p.predFilename = './pred/pishchulin13iccv/pred_sticks_lsp_oc';
        p.colorIdxs = [1 1];
    case 8
        p.name = 'Pishchulin et., ICCV''13';
        p.predFilename = './pred/pishchulin13iccv/pred_sticks_lsp_pc';
        p.colorIdxs = [1 1];
    case 9
        p.name = 'Tompson et al., NIPS''14';
        p.predFilename = './pred/tompson14nips/pred_sticks_lsp_pc';
        p.colorIdxs = [2 1];
    case 10
        p.name = 'Ramakrishna et al., ECCV''14';
        p.predFilename = './pred/ramakrishna14eccv/pred_sticks_lsp_oc';
        p.colorIdxs = [4 1];
    case 11
        p.name = 'Chen&Yuille, NIPS''14';
        p.predFilename = './pred/chen14nips/pred_sticks_lsp_oc';
        p.colorIdxs = [5 1];
    case 12
        p.name = 'Ouyang et al., CVPR''14';
        p.predFilename = './pred/ouyang14cvpr/pred_sticks_lsp_oc';
        p.colorIdxs = [7 1];
    case 13
        p.name = 'Pishchulin et al., CVPR''13';
        p.predFilename = './pred/pishchulin13cvpr/pred_sticks_lsp_oc';
        p.colorIdxs = [8 1];
    case 14
        p.name = 'Kiefel&Gehler, ECCV''14';
        p.predFilename = './pred/kiefel14eccv/pred_sticks_lsp_oc';
        p.colorIdxs = [8 2];
    case 15
        p.name = 'Kiefel&Gehler, ECCV''14';
        p.predFilename = './pred/kiefel14eccv/pred_keypoints_lsp_oc';
        p.colorIdxs = [8 3];

    case 19
        p.name = 'Bulat et al., ECCV''16';
        p.predFilename = 'pred/1st/pred_keypoints_lsp_pc';
        p.colorIdxs = [1 1];
        case 20
        p.name = 'Wei el al., CVPR''16';
        p.predFilename = 'pred/2nd/pred_keypoints_lsp_pc';
        p.colorIdxs = [2 1];
          case 21
        p.name = 'Insafutdinov et al., ECCV''16   .';
        p.predFilename = 'pred/3rd/pred_keypoints_lsp_pc';
        p.colorIdxs = [3 1];
          case 22
        p.name = 'Pishchulin et al., CVPR''16';
        p.predFilename = 'pred/4th/pred_keypoints_lsp_pc';
        p.colorIdxs = [4 1];
          case 23
        p.name = 'Lifshitz et al., ECCV''16';
        p.predFilename = 'pred/5th/pred_keypoints_lsp_pc';
        p.colorIdxs = [5 1];
          case 24
        p.name = 'Belagiannis et al., FG''17';
        p.predFilename = 'pred/6th/pred_keypoints_lsp_pc';
        p.colorIdxs = [6 1];
          case 25
        p.name = 'Yu et al., ECCV''16';
        p.predFilename = 'pred/7th/pred_keypoints_lsp_pc';
        p.colorIdxs = [7 1];
          case 26
        p.name = 'Rafi et al., BMVC''16';
        p.predFilename = 'pred/8th/pred_keypoints_lsp_pc';
        p.colorIdxs = [1 2];
          case 27
        p.name = 'Yang et al., CVPR''16';
        p.predFilename = 'pred/9th/pred_keypoints_lsp_pc';
        p.colorIdxs = [2 2];
          case 28
        p.name = 'Chen&Yuille, NIPS''14';
        p.predFilename = 'pred/10th/pred_keypoints_lsp_pc';
        p.colorIdxs = [3 2];
          case 29
        p.name = 'Fan et al., CVPR''15';
        p.predFilename = 'pred/12th/pred_keypoints_lsp_pc';
        p.colorIdxs = [4 2];
          case 30
        p.name = 'Tompson et al., NIPS''14';
        p.predFilename = 'pred/13th/pred_keypoints_lsp_pc';
        p.colorIdxs = [4 3];
          case 31
        p.name = 'Pishchulin et al., ICCV''13';
        p.predFilename = 'pred/14th/pred_keypoints_lsp_pc';
        p.colorIdxs = [4 4];
          case 32
        p.name = 'Wang&Li, CVPR''13';
        p.predFilename = 'pred/15th/pred_keypoints_lsp_pc';
        p.colorIdxs = [2 3];
       case 33  % Add GNet. Guanghan Ning. Dec 2016
        p.name = 'Ours';
        p.predFilename = 'pred/ning17iccv/pred_keypoints_lsp';
        p.colorIdxs = [8 1];

end

p.colorName = getColor(p.colorIdxs);
p.colorName = p.colorName ./ 255;

end

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose