text:
function mriCoord = cs_scs2mri(MRI, scsCoord)
% CS_SCS2MRI: Transform SCS point coordinates (in mm) to MRI coordinate system (in mm)
%
% DEPRECATED: Replace with cs_convert()

% @=============================================================================
% 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.
% =============================================================================@

disp('BST> WARNING: Deprecated function cs_scs2mri(), use cs_convert() instead.');

mriCoord = cs_convert(MRI, 'scs', 'mri', scsCoord' / 1000)' * 1000;

author: brainstorm-tools


text:
function item = cfg_exbranch(varargin)

% This is the exbranch configuration item class
%
% Data structure
% ==============
% Description fields
%    * name  - display name of config item
%    * tag   - tag of the menu item
%    * val   - 1xn cell array of cfg_items
%    * check - (optional) function handle to implement configuration
%              specific subsasgn checks based on the harvested subtree
%              rooted at this node
%    * help  - help text
% GUI/job manager fields
%    * expanded
%    * hidden
% All fields above are inherited from the branch configuration item class.
%    * prog
%    * vfiles
%    * modality
%    * vout  - function handle that generates sout struct
%    * sout  - source dependency description
%    * jout  - saved output (will be referenced by harvest of a dependency
%              target for dependency resolution at job runtime)
%    * tdeps - list where this branch is target of a dependency
%    * sdeps - list where this branch is source of a dependency
%    * chk   - field to save check status from cfg_item.check callbacks
%    * id    - id of this cfg_exbranch. This is used to reference the
%              cfg_exbranch in cfg_dep objects.
%
% Public Methods
% ==============
%    * get_strings - returns name of object
%    * gettag      - returns tag
%    * help        - returns help text
%    * harvest
%    * all_set
%
% * 'executable branch'  - See branch for details on inherited fields
%
% The layout of the configuration tree and the types of configuration items
% have been kept compatible to a configuration system and job manager
% implementation in SPM5 (Statistical Parametric Mapping, Copyright (C)
% 2005 Wellcome Department of Imaging Neuroscience). This code has been
% completely rewritten based on an object oriented model of the
% configuration tree.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_exbranch.m 5678 2013-10-11 14:58:04Z volkmar $

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

myclass = mfilename;
% Get local fields and defaults from private/mysubs_fields
[fn, defs] = mysubs_fields;

if nargin == 1
    if isstruct(varargin{1})
        % assume input is a struct to be converted back into a class
        % return with error if this does not work
        if numel(fieldnames(varargin{1})) == numel(fn)+1 && ...
                all(isfield(varargin{1}, [fn(:)' {'cfg_branch'}]))
            bitem = varargin{1}.cfg_branch;
            sitem = rmfield(varargin{1},{'cfg_branch'});
            item  = class(sitem, myclass, bitem);
            return;
        else
            cfg_message('matlabbatch:constructor:reclassify', ['Don''t know how to convert this ' ...
                                'into class ''%s''.'], myclass);
        end;
    end;
    if isa(varargin{1},myclass)
        item = varargin{1};
        return;
    end;
end;

mxpnargin = 5; % Max 5 arguments to parent initialisation
pnargin = min([nargin, mxpnargin]);
switch nargin
    case 0
        bitem = cfg_branch;
    case 1
        if isa(varargin{1},myclass)
            item = varargin{1};
            return;
        elseif isstruct(varargin{1})
            % assume input is a struct to be converted back into a class
            % no fieldname checking performed here
            bitem = varargin{1}.cfg_branch;
            sitem = rmfield(varargin{1},'cfg_branch');
            item  = class(sitem, myclass, bitem);
            return;
        else
            bitem = cfg_branch(varargin{1});
        end;
    case {2,3,4,5,6,7,8,9,10}
        bitem = cfg_branch(varargin{1:pnargin});
    otherwise
        cfg_message('matlabbatch:constructor:nargin', 'Wrong number of arguments.');
end;
for k=1:numel(fn)
    sitem.(fn{k})=defs{k};
end;
item = class(sitem, myclass, bitem);
% set additional fields (if any) - field order as in mysubs_fields
for k = 1:min(numel(fn),nargin-mxpnargin)
    item.(fn{k}) = varargin{k+mxpnargin};
end;

author: spm


text:
function geometry_test227 ( )

%*****************************************************************************80
%
%% TEST227 tests VOXELS_STEP_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST227\n' );
  fprintf ( 1, '  VOXELS_STEP_3D steps along a line from\n' );
  fprintf ( 1, '    one voxel to another.\n' );

  v1(1:dim_num) = [ 1, 1, 5 ];
  v2(1:dim_num) = v1(1:dim_num);

  inc = 7;
  jnc = 3;
  knc = -1;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  %4d  %4d  %4d  %4d\n', 0, v2(1:dim_num) );

  for i = 1 : 10
    v3 = voxels_step_3d ( v1, v2, inc, jnc, knc );
    fprintf ( 1, '  %4d  %4d  %4d  %4d\n', i, v3(1:dim_num) );
    v2(1:dim_num) = v3(1:dim_num);
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now, as a check, reverse direction and return.\n' );
  fprintf ( 1, '\n' );

  v1(1:dim_num) = v2(1:dim_num);

  inc = -inc;
  jnc = -jnc;
  knc = -knc;

  v2(1:dim_num) = v1(1:dim_num);

  fprintf ( 1, '  %4d  %4d  %4d  %4d\n', 0, v2(1:dim_num) );
  for i = 1 : 10
    v3 = voxels_step_3d ( v1, v2, inc, jnc, knc );
    fprintf ( 1, '  %4d  %4d  %4d  %4d\n', i, v3(1:dim_num) );
    v2(1:dim_num) = v3(1:dim_num);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [vertex,face] = check_face_vertex(vertex,face, options)

% check_face_vertex - check that vertices and faces have the correct size
%
%   [vertex,face] = check_face_vertex(vertex,face);
%
%   Copyright (c) 2007 Gabriel Peyre

vertex = check_size(vertex);
face = check_size(face);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a = check_size(a)
if isempty(a)
    return;
end
if size(a,1)>size(a,2)
    a = a';
end
if size(a,1)<3 && size(a,2)==3
    a = a';
end
if size(a,1)<=3 && size(a,2)>=3 && sum(abs(a(:,3)))==0
    % for flat triangles
    a = a';
end
if size(a,1)~=3 && size(a,1)~=4
    error('face or vertex is not of correct size');
end

author: gpeyre

repo: matlab-toolboxes


text:
% change in angle to closest point on wall in fly's coordinate system
function [data,units] = compute_dangle2wall(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);
  % set sign so that negative means going toward 0, positive means going
  % away from 0
  if trx(fly).nframes <= 1,
    data{i} = [];
  else
    data{i} = sign(trx(fly).angle2wall(1:end-1)).*...
      modrange(diff(trx(fly).angle2wall,1,2),-pi,pi)./trx(fly).dt;
  end
end
units = parseunits('rad/s');

author: kristinbranson

repo: JAABA


text:
function varargout = resolutionmerge(varargin)
% RESOLUTIONMERGE M-file for resolutionmerge.fig
%
% GUI for improving resolution of lower resolution image using higher
% resolution image using RGB-to-HSI conversion.
% Images have to be spatially registered.
%
% Monochromatic (gray scale) low resolution image is converged to colored
% image by false color mapping to "hot" color scheme. This RGB image is
% then converted to Hue, Saturation and Value (HSV) image. Value component is replaced
% by higher resolution image and the resulting HSV image is converted back to
% RGB. This RGB merged image converted to gray scale is a merged image with
% improved spatial resolution.
%
% To run, type:
% >> resolutionmerge
%
% Can work with any type of color or monochromatic image. If images are colored
% they are converted to gray scale first and then merged.
% To load variables from matlab file, the file should contain variables
% LOWRES and HIGHRES
%
% created by K.Artyushkova
% February 2010
%
%

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @resolutionmerge_OpeningFcn, ...
                   'gui_OutputFcn',  @resolutionmerge_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 resolutionmerge is made visible.
function resolutionmerge_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 resolutionmerge (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

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

% --- Outputs from this function are returned to the command line.
function varargout = resolutionmerge_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;

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

[filename, pathname]=uigetfile('*.*','Open low resolution image');
cd(pathname)
[N,M]=size(filename);
image=imread(char(filename));
[n,m,p]=size(image);
if p==3
    lowres=rgb2gray(image);
else
    lowres=image;
end
axes(handles.axes1)
handles.lowres=lowres;
iptsetpref('ImshowAxesVisible', 'on')
imagesc(uint8(lowres)), colormap(gray)
guidata(hObject,handles)

% --------------------------------------------------------------------
function open_lowres_Callback(hObject, eventdata, handles)
% hObject    handle to open_lowres (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename, pathname]=uigetfile('*.mat','Open low resolution image. Image should be saved in variable LOWRES');
cd(pathname)
D=load (filename, 'lowres');
image=D.lowres;
[n,m,p]=size(image);
if p==3
    lowres=rgb2gray(image);
else
    lowres=image;
end
axes(handles.axes1)
handles.lowres=lowres;
iptsetpref('ImshowAxesVisible', 'on')
imagesc(uint8(lowres)), colormap(gray)
guidata(hObject,handles)

% --------------------------------------------------------------------
function load_highres_Callback(hObject, eventdata, handles)
% hObject    handle to load_highres (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename, pathname]=uigetfile('*.*','Open high resolution image');
cd(pathname)
[N,M]=size(filename);
image=imread(char(filename));
[n,m,p]=size(image);
if p==3
    highres=rgb2gray(image);
else
    highres=image;
end
axes(handles.axes2)
handles.highres=highres;
iptsetpref('ImshowAxesVisible', 'on')
imagesc(uint8(highres)), colormap(gray)
guidata(hObject,handles)

% --------------------------------------------------------------------
function Open_highres_Callback(hObject, eventdata, handles)
% hObject    handle to Open_highres (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename, pathname]=uigetfile('*.mat','Open high resolution image. Image should be saved in variable LOWRES');
cd(pathname)
D=load (filename, 'highres');
image=D.highres;
[n,m,p]=size(image);
if p==3
    highres=rgb2gray(image);
else
    highres=image;
end
axes(handles.axes2)
handles.highres=highres;
iptsetpref('ImshowAxesVisible', 'on')
imagesc(uint8(highres)), colormap(gray)
guidata(hObject,handles)

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

%% Main algorithm
% --- Executes on button press in Merge_resolution.
function Merge_resolution_Callback(hObject, eventdata, handles)
% hObject    handle to Merge_resolution (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
lowres=handles.lowres;
highres=handles.highres;
map=hot(256);
R=ind2rgb(lowres,map);
axes(handles.axes3)
iptsetpref('ImshowAxesVisible', 'on')
imagesc(R), colormap(hot)
HSV=rgb2hsv(R);
HSV_sub=HSV;
HSV_sub(:,:,3)=highres;
merged=hsv2rgb(HSV_sub);
merged=uint8(merged);
mergedG=rgb2gray(merged);
axes(handles.axes4)
iptsetpref('ImshowAxesVisible', 'on')
imagesc(uint8(mergedG)), colormap(gray)
handles.merged=merged;
handles.mergedG=mergedG;
guidata(hObject, handles);

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

% Hint: get(hObject,'Value') returns toggle state of displaycolor
merged=handles.merged;
axes(handles.axes4)
iptsetpref('ImshowAxesVisible', 'on')
imagesc(uint8(merged))

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

% Hint: get(hObject,'Value') returns toggle state of disp_grayscale
mergedG=handles.mergedG;
axes(handles.axes4)
iptsetpref('ImshowAxesVisible', 'on')
imagesc(uint8(mergedG)), colormap(gray)

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

%% Output section
% --------------------------------------------------------------------
function save_mat_Callback(hObject, eventdata, handles)
% hObject    handle to save_mat (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
mergedG=handles.mergedG;
merged=handles.merged;
datapath = uigetdir;
cd(datapath)
[filename, pathname] = uiputfile('*.mat', 'save merged image into mat file');
save(filename, 'merged','mergedG')

% --------------------------------------------------------------------
function save_tiff_Callback(hObject, eventdata, handles)
% hObject    handle to save_tiff (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
mergedG=handles.mergedG;
merged=handles.merged;
datapath = uigetdir;
cd(datapath)
imwrite(uint8(mergedG),'mergedG.tiff','tiff')
imwrite(uint8(merged),'merged.tiff','tiff')

author: Sable

repo: mcbench-benchmarks


text:
function [ patches, pdm, clmParams ] = Load_CLNF_general()
%LOAD_CLNF_GENERAL Summary of this function goes here
%   Detailed explanation goes here
    clmParams = struct;

    clmParams.window_size = [25,25; 23,23; 21,21];
    clmParams.numPatchIters = size(clmParams.window_size,1);

    [patches] = Load_Patch_Experts( '../models/general/', 'ccnf_patches_*_general.mat', [], [], clmParams);

    % the default PDM to use
    pdmLoc = ['../models/pdm/pdm_68_aligned_wild.mat'];
    load(pdmLoc);

    pdm = struct;
    pdm.M = double(M);
    pdm.E = double(E);
    pdm.V = double(V);

    % the default model parameters to use
    clmParams.regFactor = [35, 27, 20];
    clmParams.sigmaMeanShift = [1.25, 1.375, 1.5];
    clmParams.tikhonov_factor = [2.5, 5, 7.5];

    clmParams.startScale = 1;
    clmParams.num_RLMS_iter = 10;
    clmParams.fTol = 0.01;
    clmParams.useMultiScale = true;
    clmParams.use_multi_modal = 1;
    clmParams.multi_modal_types  = patches(1).multi_modal_types;

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function D = spm_eeg_inv_fmripriors(S)
% Generate fMRI priors for the M/EEG source reconstruction
% FORMAT D = spm_eeg_inv_fmripriors(S)
%
% S        - optional input struct
% (optional) fields of S:
%  .D      - MEEG object or filename of M/EEG mat-file
%  .fmri   - filename of prior (SPM) image to be used
%  [.gm    - filename of grey matter (GM) image] {unused}
%  .space  - native (0) or MNI (1) space (must be same for SPM and GM images)
%  .hthr   - height threshold of prior image [defaults: 0.5]
%  .ethr   - extent threshold of clusters in prior image [default: 1]
%  .ncomp  - maximal number of priors component to be extracted [default: Inf]
%  .smooth - variance of the smoothing kernel onto the surface [default: 0.2] {unused}
%  .disp   - whether to display priors on mesh [default: 0]
%
% D.inv{D.val}.inverse.fmri.priors   - MAT filename containing a variable 'pQ' that
%            is a [ncomp] cell array of [nb vertices] vectors describing spatial priors
% D.inv{D.val}.inverse.fmri.texture  - GIfTI texture filename containing all
%            spatial priors
% D.inv{D.val}.inverse.fmri.clusters - image filename containing clusters as labels
%__________________________________________________________________________
%
% Reference:
%
% A Parametric Empirical Bayesian framework for fMRI-constrained MEG/EEG
% source reconstruction. Henson R, Flandin G, Friston K & Mattout J.
% Human Brain Mapping (in press).
%__________________________________________________________________________
% Copyright (C) 2008-2015 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin and Rik Henson
% $Id: spm_eeg_inv_fmripriors.m 6301 2015-01-12 17:23:08Z guillaume $

%-Get MEEG object
%--------------------------------------------------------------------------
try
    D = S.D;
catch
    [D, sts]   = spm_select(1, 'mat', 'Select M/EEG mat file');
    if ~sts, D = []; return; end
    S.D = D;
end

[D, val] = spm_eeg_inv_check(D);

%-Input parameters
%--------------------------------------------------------------------------
try
    S.fmri;
catch
    [S.fmri, sts] = spm_select(1,'image','Select prior (eg thresholded SPM) image...');
    if ~sts, return; end
end

% try
%     S.gm;
% catch
%     [S.gm, sts] = spm_select([0 1],'image','Select grey matter (GM) image...');
% end

try
    space = S.space;   % 0 = native, 1 = MNI
catch
    space = spm_input('Image space','+1','b',{'Native|MNI'},[0 1],1);
end

if space
    try
        m = D.inv{val}.mesh.tess_mni;
        %fprintf('Note: assuming SPM and GM images are in MNI space...\n');
    catch
        error('This D structure has no MNI cortical mesh stored.');
    end
else
    try
        m = export(gifti(D.inv{val}.mesh.tess_ctx),'spm');
        %fprintf('Note: assuming SPM and GM images are in subject''s native MRI space...\n');
    catch
        error('This D structure has no cortical mesh stored.');
    end
end

try
    S.hthr;
catch
    S.hthr = 0.5; % assume binary prior
end

try
    S.ethr;
catch
    S.ethr = 1; % no threshold on extent
end

try
    S.ncomp;
catch
    S.ncomp = Inf; % all components
end

try
    S.bincomp;
catch
    S.bincomp = 1; % default to binary priors
end

try
    S.varcomp;
catch
    S.varcomp = 1; % default to variance priors (vectors)
end

% try
%     S.smooth;
% catch
%     S.smooth = 0.2;
% end

try
    S.disp;
catch
    S.disp = 0;
end

%-Extracting clusters from functional image
%==========================================================================
V     = spm_vol(S.fmri);
prior = spm_read_vols(V);

%-Height threshold
%--------------------------------------------------------------------------
prior = prior > S.hthr;

%-Connected Component labelling
%--------------------------------------------------------------------------
[l2, num] = spm_bwlabel(double(prior),26);
if ~num
    fprintf('No suprathreshold clusters available.\n');
    return
end

%-Extent threshold, and sort clusters according to their extent
%--------------------------------------------------------------------------
[n, ni] = sort(histc(l2(:),0:num), 1, 'descend');
l  = zeros(size(l2));
n  = n(2:end);      ni = ni(2:end)-1;
ni = ni(n>=S.ethr); n  = n(n>=S.ethr);
S.ncomp = min(S.ncomp, length(n));
for i=1:S.ncomp
    l(l2==ni(i)) = i;
end
clear l2 ni
fprintf('Selected %d clusters (out of %d) in prior image.\n',S.ncomp,num);

%-Projecting volumetric clusters on surface mesh
%==========================================================================
q = zeros(S.ncomp, size(m.vert,1));
for i=1:S.ncomp
    q(i,:) = spm_mesh_project(m.vert,struct('dat',double(l==i),'mat',V.mat),'nn');
end
q(~any(q,2),:) = [];
fprintf('After projection, %d clusters remaining.\n',size(q,1));
if isempty(q), return; end

%-Smooth, binarize and save in output variable
%--------------------------------------------------------------------------
pQ = cell(1,size(q,1));
for i = 1:size(q,1)
    qq    = q(i,:)';
    %qq   = spm_mesh_smooth(struct('faces',double(m.face),'vertices',m.vert),qq, S.smooth);
    qq    = qq .* (qq > exp(-8));
    if S.bincomp
        pQ{i} =  double(qq > 0)';  % binarise
    end
end

%-Display and export clusters
%==========================================================================
if S.disp
    for i=1:numel(pQ)
        spm_eeg_render(struct('faces',double(m.face),'vertices',m.vert),...
            struct('texture',pQ{i}));
    end
end

%-Save clusters as an image of labels
%--------------------------------------------------------------------------
[pth,name] = fileparts(S.fmri);
D.inv{val}.inverse.fmri.clusters = fullfile(pth,['cluster_' name spm_file_ext]);
V = struct('fname',   D.inv{val}.inverse.fmri.clusters, ...
           'dim',     V.dim, ...
           'dt',      [spm_type('uint16') spm_platform('bigend')], ...
           'mat',     V.mat, ...
           'pinfo',   [1 0 0]', ...
           'descrip', 'clusters');
V = spm_write_vol(V,l);

%-Save spatial priors vectors as GIfTI file
%--------------------------------------------------------------------------
[pth,name] = fileparts(D.fname);
D.inv{val}.inverse.fmri.texture = fullfile(pth,['priors_' name '_' num2str(val) '.func.gii']);
G          = gifti;
G.cdata    = cat(1, pQ{:})';
save(G,D.inv{val}.inverse.fmri.texture);

%-Save spatial priors vectors as MAT-file
%--------------------------------------------------------------------------
[pth,name] = fileparts(D.fname);
D.inv{val}.inverse.fmri.priors = fullfile(pth,['priors_' name '_' num2str(val) '.mat']);
if ~S.varcomp
    pQ = pQ*pQ';
end

save(D.inv{val}.inverse.fmri.priors,'pQ', spm_get_defaults('mat.format'));

%-Save D structure
%--------------------------------------------------------------------------
%D.save;


text:
%% Example
%
% Copyright:
%   2018 - Simo Särkkä and Arno Solin
%
% License:
%   This software is provided under the MIT License. See the accompanying
%   LICENSE file for details.

%% Leapfrog / Verlet integration

dts = logspace(-2,0,10);
err = nan(10,numel(dts));

for j=1:numel(dts)
    fprintf('Running step %d/%d\n',j,numel(dts));

% Leapfrog / Verlet integration

  % Lock seed
  if exist('rng') % Octave doesn't have rng
    rng(j,'twister')
  else
    randn('state',j);
  end

  % Parameters
  n = 25000;
  g = 1; % g>0
  eta = 1;
  q = 1;
  x0 = 1;
  v0 = 0;

  % Time discretization
  dt = dts(j); %0.1;
  t = 0:dt:10;

  % Specify the model
  f = @(x) -g*x;
  s = @(x) 1;

  % Allocate space and set initial conditions
  x = zeros(n,numel(t)); x(:,1) = x0;
  v = zeros(n,numel(t)); v(:,1) = v0;

  % For each trajectory
  for i=1:n

    % Simulate one trajectory from the leapfrog method
    z = leapfrog(f,s,eta,t,[x0; v0],q);

    x(i,:) = z(1,:);
    v(i,:) = z(2,:);

  end

  figure(1); clf
  subplot(211)
    plot(t,x(1:5,:), ...
         t,mean(x), ...
         t,mean(x)+1.96*std(x),'--', ...
         t,mean(x)-1.96*std(x),'--')
    xlabel('Time, t'), ylabel('Position, x(t)')
  subplot(212)
    plot(t,v(1:5,:), ...
         t,mean(v), ...
         t,mean(v)+1.96*std(v),'--', ...
         t,mean(v)-1.96*std(v),'--')
    xlabel('Time, t'), ylabel('Velocity, v(t)')

% Euler-Maruyama

  % Same random seed
  % Lock seed
  if exist('rng') % Octave doesn't have rng
    rng(j,'twister')
  else
    randn('state',j);
  end

  % Specify the model
  f = @(x,t) [0 1; -g -eta]*x;
  L = @(x,t) [0; 1];

  % Allocate space and set initial conditions
  x_em = zeros(n,numel(t)); x_em(:,1) = x0;
  v_em = zeros(n,numel(t)); v_em(:,1) = v0;

  % For each trajectory
  for i=1:n

    % Simulate one trajectory from the leapfrog method
    z = eulermaruyama(f,L,t,[x0; v0],q);

    x_em(i,:) = z(1,:);
    v_em(i,:) = z(2,:);

  end

  figure(1); clf
  subplot(211)
    plot(t,x_em(1:5,:), ...
         t,mean(x_em), ...
         t,mean(x_em)+1.96*std(x_em),'--', ...
         t,mean(x_em)-1.96*std(x_em),'--')
    xlabel('Time, t'), ylabel('Position, x(t)')
  subplot(212)
    plot(t,v(1:5,:), ...
         t,mean(v_em), ...
         t,mean(v_em)+1.96*std(v_em),'--', ...
         t,mean(v_em)-1.96*std(v_em),'--')
    xlabel('Time, t'), ylabel('Velocity, v(t)')

% Closed-form solution

  % This is a LTI SDE of the form
  F = [0 1; -g -eta];
  L = [0; 1];
  Qc = q;

  % Discretize
  [A,Q] = lti_disc(F,L,Qc,dt);

  % Evaluate
  M = zeros(2,numel(t)); M(:,1) = [x0; v0];
  P = zeros(2,2,numel(t));

  % Loop through
  for k=1:numel(t)-1
    M(:,k+1) = A*M(:,k);
    P(:,:,k+1) = A*P(:,:,k)*A'+Q;
  end

  figure(2); clf
  subplot(211)
    plot(t,M(1,:), ...
         t,M(1,:)'+1.96*sqrt(squeeze(P(1,1,:))),'--', ...
         t,M(1,:)'-1.96*sqrt(squeeze(P(1,1,:))),'--')
  subplot(212)
    plot(t,M(2,:), ...
         t,M(2,:)'+1.96*sqrt(squeeze(P(2,2,:))),'--', ...
         t,M(2,:)'-1.96*sqrt(squeeze(P(2,2,:))),'--')

% Calculate errors

  err(1,j) = P(1,1,end)-std(x(:,end))^2;
  err(2,j) = P(2,2,end)-std(v(:,end))^2;

  err(3,j) = P(1,1,end)-std(x_em(:,end))^2;
  err(4,j) = P(2,2,end)-std(v_em(:,end))^2;

  err(5,j) = M(1,end)-mean(x(:,end));
  err(6,j) = M(2,end)-mean(v(:,end));

  err(7,j) = M(1,end)-mean(x_em(:,end));
  err(8,j) = M(2,end)-mean(v_em(:,end));

  err(9,j)  = mean(abs(squeeze(P(1,1,:))'-std(x).^2));
  err(10,j) = mean(abs(squeeze(P(1,1,:))'-std(x_em).^2));

end

%% Compose results

  figure(1); clf
    loglog(dts,abs(err(9,:)),'-k', ...
           dts,abs(err(10,:)),'--k')
    xlabel('Time step length, $\Delta t$')
    ylabel('Mean absolute error in $\sigma_x^2$')
    legend('Leapfrog Verlet', 'Euler--Maruyama')
    %ylim([0.4 10])
    xlim([min(dts) max(dts)])
    set(gca,'XTick',[0.01 0.1 1],'XTickLabel',[0.01 0.1 1])
    set(gca,'YTick',[1e-3 1e-2 1e-1 1e0 1e1], ...
            'YTickLabel',{'','$10^{-2}$','$10^{-1}$','$10^{0}$','$10^1$'})

author: AaltoML

repo: SDE

save_path: github-repos/MATLAB/AaltoML-SDE


text:
function [ errors ] = test_operators( )

errors = 0;

errors = errors + test_fourier_basis_ring();
errors = errors + test_fourier_basis_ring_advanced();

%errors = errors + test_translate();
errors = errors + test_modulate();

errors = errors + test_gwft();
errors = errors + test_gft();
errors = errors + test_ngwft();

errors = errors + test_divgradient();
errors = errors + test_gradient();
errors = errors + test_gradient2();
errors = errors + test_gradient_mat();

errors = errors + test_directed();
errors = errors + test_advection();

errors = errors + test_div();
errors = errors + test_div2();
errors = errors + test_div3();
errors = errors + test_laplacian();

errors = errors + test_kron_reduction();

end

function s = eigvec_ring(input,N)
s = zeros(N,numel(input));
n = (0:(N-1))';
for ind = 1:numel(input)
    if mod(input(ind),2) % odd
        s(:,ind) = 1/sqrt(N)*exp(2*pi*1i*(N-(input(ind)+1)/2)*n/N);
    else % even
        s(:,ind) = 1/sqrt(N)*exp(pi*1i*input(ind)*n/N);
    end
end

end

function s = eigval_ring(input,N)
s = zeros(size(input));
for ind = 1:numel(input)
    if mod(input(ind),2) % odd
        s(ind) = 1-cos(pi * (input(ind)+1)/N );
    else % even
        s(ind) = 1-cos(pi * (input(ind))/N );
    end
end

end

function errors = test_fourier_basis_ring_advanced()
errors = 0;
N = 10;
e = 2*eigval_ring(0:N-1,N)';
U = eigvec_ring(0:N-1,N);
G = gsp_ring(N);
G = gsp_compute_fourier_basis(G);
errors = errors +gsp_assert_test(G.e,e,1e-14,'test ring eigenvalues N=10');
errors = errors +gsp_assert_test(G.U,U,1e-14,'test ring eigenvector N=10');

N = 11;
e = 2*eigval_ring(0:N-1,N)';
U = eigvec_ring(0:N-1,N);
G = gsp_ring(N);
G = gsp_compute_fourier_basis(G);
errors= errors +gsp_assert_test(G.e,e,1e-14,'test ring eigenvalues N=11');
errors = errors +gsp_assert_test(G.U,U,1e-14,'test ring eigenvector N=11');

end

function errors = test_laplacian()
errors = 0;
N = 10;
G = gsp_path(N);

x = rand(N,1);

y1 = G.L*x;
y2 = -laplacianx_op(x);

if norm(y1-y2)/norm(y1)>10e-10;
   errors = errors +1;
   norm(y1-y2)/norm(y1)
   warning('OPERATOR: Error in laplacian test 1')

else
    fprintf('OPERATOR: laplacian 1 ok\n');
end

end

function errors = test_div()
errors = 0;
N = 10;
G = gsp_path(N);
G = gsp_adj2vec(G);

x = rand(N-1,1);

y1 = gsp_div(G,x);
y2 = -div_op1d([x;0]);

if norm(y1-y2)/norm(y1)>10e-10;
   errors = errors +1;
   norm(y1-y2)/norm(y1)
   warning('OPERATOR: Error in div test 1')

else
    fprintf('OPERATOR: div 1 ok\n');
end

end

function errors = test_div2()
errors = 0;
N = 10;
G = gsp_sensor(N);
G = gsp_adj2vec(G);

M = 5;

x = rand(N,M);

t = gsp_grad(G,x);

y1 = gsp_div(G,t);
y2 = gsp_div_old(G,t);

if norm(y1(:)-y2(:))/norm(y1(:))>10e-10;
   errors = errors +1;
   norm(y1(:)-y2(:))/norm(y1(:))
   warning('OPERATOR: Error in div test 2')

else
    fprintf('OPERATOR: div 2 ok\n');
end

end

function errors = test_div3()
errors = 0;
N = 10;
G = gsp_sensor(N);
G = gsp_adj2vec(G);

M = 5;

x = rand(N,M);

t = gsp_grad(G,x);

y1 = gsp_div(G,t);

y2 = zeros(N,M);
for ii = 1:M
    y2(:,ii) = gsp_div_old(G,t(:,ii));
end

if norm(y1(:)-y2(:))/norm(y1(:))>10e-10;
   errors = errors +1;
   norm(y1(:)-y2(:))/norm(y1(:))
   warning('OPERATOR: Error in div test 3')

else
    fprintf('OPERATOR: div 3 ok\n');
end

end

function errors = test_gradient()
errors = 0;
N = 50;
G = gsp_path(N);
G = gsp_adj2vec(G);

x = rand(N,1);

y1 = gsp_grad(G,x);
y2 = gradient_op1d(x);

if norm(y1-y2(1:(end-1)))/norm(y1)>10e-10;
   errors = errors +1;
   norm(y1-y2(1:(end-1)))/norm(y1)
   warning('OPERATOR: Error in gradient test 1')

else
    fprintf('OPERATOR: gradient 1 ok\n');
end

end

function errors = test_gradient2()
errors = 0;
N = 50;
G = gsp_sensor(N);
G = gsp_adj2vec(G);

M = 5;

x = rand(N,M);

y1 = gsp_grad(G,x);

y2 = zeros(length(G.v_in),M);
for ii = 1:M
    y2(:,ii) = gsp_grad(G,x(:,ii));
end

if norm(y1(:)-y2(:))/norm(y1(:))>10e-10;
   errors = errors +1;
   norm(y1(:)-y2(:))/norm(y1(:))
   warning('OPERATOR: Error in gradient test 2')

else
    fprintf('OPERATOR: gradient 2 ok\n');
end

end

function errors = test_gradient_mat()
errors = 0;
N = 50;
G = gsp_sensor(N);
G = gsp_adj2vec(G);

x = rand(N,1);

y1 = gsp_grad(G,x);

D = gsp_grad_mat(G);

y2 = D*x;

if norm(y1(:)-y2(:))/norm(y1(:))>10e-10;
   errors = errors +1;
   norm(y1(:)-y2(:))/norm(y1(:))
   warning('OPERATOR: Error in gradient mat test')

else
    fprintf('OPERATOR: gradient mat ok\n');
end

end

function errors = test_directed()
errors = 0;

try

    G = gsp_sensor(50);

    Gd = gsp_assign_rand_direction(G,.5);
    Gd = gsp_adj2vec(Gd);

    fprintf('OPERATOR: directed ok\n');
catch
    errors = errors + 1;
    warning('OPERATOR: Error directed test')
end

end

function errors = test_advection()
errors = 0;
G = gsp_sensor();

Gd = gsp_assign_rand_direction(G,.5);
Gd = gsp_adj2vec(Gd);

W = diag(diag(Gd.Adv)) - Gd.Adv';

err = norm(Gd.W-W,'fro')/norm(Gd.W,'fro');

if err>1e-10;
   errors = errors +1;
   err
   warning('OPERATOR: Error advection test')

else
    fprintf('OPERATOR: advection ok\n');
end

end

function errors = test_divgradient()
errors = 0;
N = 50;
G = gsp_sensor(N);
G = gsp_adj2vec(G);

x = rand(N,1);

y1 = gsp_div(G,gsp_grad(G,x));
y2 = G.L*x;

if norm(y1-y2)/norm(y1)>10e-10
   errors = errors +1;
   norm(y1-y2)/norm(y1)
   warning('OPERATOR: Error in divgradient test 1')

else
    fprintf('OPERATOR: divgradient 1 ok\n');
end

if norm(full(G.Diff'*G.Diff-G.L))>10e-10;
   errors = errors +1;
   norm(full(G.Diff'*G.Diff-G.L))
   warning('OPERATOR: Error in divgradient test 2')

else
    fprintf('OPERATOR: divgradient 2 ok\n');
end

G = gsp_create_laplacian(G,'normalized');
G = gsp_adj2vec(G);

y1 = gsp_div(G,gsp_grad(G,x));
y2 = G.L*x;

if norm(y1-y2)/norm(y1)>10e-10
   errors = errors +1;
   norm(y1-y2)/norm(y1)
   warning('OPERATOR: Error in divgradient test 3')

else
    fprintf('OPERATOR: divgradient 3 ok\n');
end

if norm(full(G.Diff'*G.Diff-G.L))>10e-10;
   errors = errors +1;
   norm(full(G.Diff'*G.Diff-G.L))
   warning('OPERATOR: Error in divgradient test 3')

else
    fprintf('OPERATOR: divgradient 3 ok\n');
end

end

function errors = test_fourier_basis_ring()
errors = 0;
N = 50;
G = gsp_ring(N);
G = gsp_compute_fourier_basis(G);

if sum(diff(G.e) < -1e-10)
   errors = errors +1;
   warning('OPERATOR: Error in the Fourier test 1')

else
    fprintf('OPERATOR: Fourier 1 ok\n');
end

if sum(norm(G.L-G.U*diag(G.e)*G.U','fro') < -1e-10)
   errors = errors +1;
   warning('OPERATOR: Error in the Fourier test 2')

else
    fprintf('OPERATOR: Fourier 2 ok\n');
end

N = 51;
G = gsp_ring(N);
G = gsp_compute_fourier_basis(G);

if sum(norm(G.L-G.U*diag(G.e)*G.U','fro') < -1e-10)
   errors = errors +1;
   warning('OPERATOR: Error in the Fourier test 2')

else
    fprintf('OPERATOR: Fourier 2 ok\n');
end

end

function errors = test_translate()
errors = 0;
N = 50;
k = 10;
G = gsp_ring(N);
G = gsp_compute_fourier_basis(G);

f = rand(N,1);

ft = gsp_translate(G,f,k);

if norm(ft-circshift(f,[-k+1,0]))<1e-10
    fprintf('OPERATOR: translate ok\n');
else
    errors = errors +1;
   warning('OPERATOR: Error in the translate test')

end

end

function errors = test_modulate()
errors = 0;
N = 30;
k = 10;
G = gsp_ring(N);
G = gsp_compute_fourier_basis(G);

f = rand(N,1);

% compute the ordering vectors
v = gsp_classic2graph_eig_order( N );
iv = 1:N;
iv(v) = iv;

ft = gsp_modulate(G,f,iv(k));

fthat(v) = gsp_gft(G,ft);
ft2hat(v) = gsp_gft(G,f); % back into the classical ordering
ft2hat = circshift(ft2hat,[0,-k+1]); % circshift

if norm(fthat-ft2hat)<1e-10
    fprintf('OPERATOR: modulate ok\n');
else
    errors = errors +1;
   warning('OPERATOR: Error in the modulate test')

end

end

function errors = test_kron_reduction()

errors = 0;

try
    figure(100)
      N = 64;
      param.distribute = 1;
      param.Nc = 5;
      param.regular = 1;
      G = gsp_sensor(N,param);
      ind = 1:2:N;
      Gnew = gsp_kron_reduce( G,ind );
      subplot(121)
      gsp_plot_graph(G);
      title('Original graph');
      subplot(122)
      gsp_plot_graph(Gnew);
      title('Kron reduction');

    close(100)

    fprintf('OPERATOR: KRON ok\n');
catch
    errors = errors +1;
   warning('OPERATOR: Error in the KRON test')

end

gsp_reset_seed(1)
      G1 = gsp_kron_reduce( G,ind );
      gsp_reset_seed(1)
      Lnew= gsp_kron_reduce_old( G.L,ind );
if norm(Lnew-G1.L,'fro')<1e-10
      fprintf('OPERATOR: KRON ok\n');
else
          errors = errors +1;
    norm(G2.W-G1.W,'fro')
   warning('OPERATOR: Error in the KRON test')
end

end

function errors = test_gft()

N = 30;
G = gsp_sensor(N);
G = gsp_compute_fourier_basis(G);

X = rand(N,20,10);

Xhat1 = gsp_gft(G,X);

Xhat2 = zeros(N,20,10);
for ii = 1:10
   Xhat2(:,:,ii) = gsp_gft(G, X(:,:,ii));
end

errors = gsp_assert_test(Xhat1,Xhat2,1e-10,'GFT');

end

repo: gspbox


text:
function imgHiRes = sr_demo(filePath, opt)
% SR_DEMO
%
% Input:
%     - filePath: filePath contains
%       (1) dataPath: path to the input images
%       (2) imgFileName: the file name of the input low-resolution image
%     - opt: algorithm parameters, obtained by
%            opt = sr_init_opt(SRF);
% Output:
%     - imgHiRes: the super-resolved images
%
% Example usage:
%   filePath.dataPath    = 'data\Urban100\image_SRF_4';
%   filePath.imgFileName = 'img_001_SRF_4_LR.png';
%
%   imgHiRes = sr_demo(filePath, opt);
% Disclaimer:
%
%   It is provided for educational/researrch purpose only.
%   If you find the software useful, please consider cite our paper.
%
%   Single Image Super-Resolution using Transformed Self-Exemplars
%   Jia-Bin Huang, Abhishek Singh, and Narendra Ahuja
%   IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2015
%
% Contact:
%   Jia-Bin Huang
%   Electrical and Computer Engineering
%   University of Illinois, Urbana-Champaign
%   www.jiabinhuang.com

% =========================================================================
% Extract planar structures
% =========================================================================
fprintf('- Extract planar structures \n');
tic;
modelPlane = sr_extract_plane(filePath.dataPath, filePath.imgFileName, opt);
tAnalysis = toc;
fprintf('Done in %6.3f seconds.\n\n', tAnalysis);

% =========================================================================
% Construct image pyramid for coarse-to-fine image super-resolution
% =========================================================================
tic;
% Load the input low-resolution image:
img = im2single(imread(fullfile(filePath.dataPath, filePath.imgFileName)));
if(ndims(img)~=3) % support only RGB images
    img = cat(3, img, img, img);
end
tImgPyramid = toc;
fprintf('Done in %6.3f seconds.\n\n', tImgPyramid);

% Create image pyramid and multi-level planar structure constraints
fprintf('- Construct plane pyramid: \n');
tic;
[imgPyrH, imgPyrL, scaleImgPyr] = sr_create_img_pyramid(img, opt);
modelPlane = sr_planar_structure_pyramid(scaleImgPyr, modelPlane, opt.topLevel);
tPlanePyramid = toc;
fprintf('Done in %6.3f seconds.\n\n', tPlanePyramid);

% =========================================================================
% Super-resolution by patch-based synthesis
% =========================================================================
fprintf('- Single Image Super-Resolution using Transformed Self-Exemplars \n');
tic;
imgPyr = sr_synthesis(imgPyrH, imgPyrL, scaleImgPyr, modelPlane, filePath, opt);
tSynthesis = toc;
fprintf('Synthesis took %6.3f seconds.\n', tSynthesis);

% Get the level corresponding to the desired high-resolution image
if(mod(opt.SRF,2) == 0 )
    lvlInd = opt.origResLvl - log(opt.SRF)/log(2)*opt.nLvlToRedRes;
else
    lvlInd = opt.origResLvl - log(opt.SRF)/log(3)*opt.nLvlToRedRes;
end

imgHiRes = imgPyr{uint8(lvlInd)};

% Visualize plane parameters
% if(0)
%     imgTformVis = sr_vis_planar_tform(imgPyr{opt.origResLvl-6}, modelPlane{opt.origResLvl-6});
%
%     imgVisResName = [filePath.imgFileName(1:end-4), '_vis.png'];
%     imwrite(imgTformVis, fullfile('result\vis', imgVisResName));
% end

end

repo: SelfExSR


text:
function [mssim, ssim_map] = ssim_index(img1, img2, K, window, L)

%========================================================================
%SSIM Index, Version 1.0
%Copyright(c) 2003 Zhou Wang
%All Rights Reserved.
%
%The author was with Howard Hughes Medical Institute, and Laboratory
%for Computational Vision at Center for Neural Science and Courant
%Institute of Mathematical Sciences, New York University, USA. He is
%currently with Department of Electrical and Computer Engineering,
%University of Waterloo, Canada.
%
%----------------------------------------------------------------------
%Permission to use, copy, or modify this software and its documentation
%for educational and research purposes only and without fee is hereby
%granted, provided that this copyright notice and the original authors'
%names appear on all copies and supporting documentation. This program
%shall not be used, rewritten, or adapted as the basis of a commercial
%software or hardware product without first obtaining permission of the
%authors. The authors make no representations about the suitability of
%this software for any purpose. It is provided "as is" without express
%or implied warranty.
%----------------------------------------------------------------------
%
%This is an implementation of the algorithm for calculating the
%Structural SIMilarity (SSIM) index between two images. Please refer
%to the following paper:
%
%Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image
%quality assessment: From error measurement to structural similarity"
%IEEE Transactios on Image Processing, vol. 13, no. 4, Apr. 2004.
%
%Kindly report any suggestions or corrections to zhouwang@ieee.org
%
%----------------------------------------------------------------------
%
%Input : (1) img1: the first image being compared
%        (2) img2: the second image being compared
%        (3) K: constants in the SSIM index formula (see the above
%            reference). defualt value: K = [0.01 0.03]
%        (4) window: local window for statistics (see the above
%            reference). default widnow is Gaussian given by
%            window = fspecial('gaussian', 11, 1.5);
%        (5) L: dynamic range of the images. default: L = 255
%
%Output: (1) mssim: the mean SSIM index value between 2 images.
%            If one of the images being compared is regarded as
%            perfect quality, then mssim can be considered as the
%            quality measure of the other image.
%            If img1 = img2, then mssim = 1.
%        (2) ssim_map: the SSIM index map of the test image. The map
%            has a smaller size than the input images. The actual size:
%            size(img1) - size(window) + 1.
%
%Default Usage:
%   Given 2 test images img1 and img2, whose dynamic range is 0-255
%
%   [mssim ssim_map] = ssim_index(img1, img2);
%
%Advanced Usage:
%   User defined parameters. For example
%
%   K = [0.05 0.05];
%   window = ones(8);
%   L = 100;
%   [mssim ssim_map] = ssim_index(img1, img2, K, window, L);
%
%See the results:
%
%   mssim                        %Gives the mssim value
%   imshow(max(0, ssim_map).^4)  %Shows the SSIM index map
%
%========================================================================

if (nargin < 2 | nargin > 5)
   mssim = -Inf;
   ssim_map = -Inf;
   return;
end

if (size(img1) ~= size(img2))
   mssim = -Inf;
   ssim_map = -Inf;
   return;
end

[M N] = size(img1);

if (nargin == 2)
   if ((M < 11) | (N < 11))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   window = fspecial('gaussian', 11, 1.5);	%
   K(1) = 0.01;								      % default settings
   K(2) = 0.03;								      %
   L = 255;                                  %
end

if (nargin == 3)
   if ((M < 11) | (N < 11))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   window = fspecial('gaussian', 11, 1.5);
   L = 255;
   if (length(K) == 2)
      if (K(1) < 0 | K(2) < 0)
		   mssim = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   mssim = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

if (nargin == 4)
   [H W] = size(window);
   if ((H*W) < 4 | (H > M) | (W > N))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   L = 255;
   if (length(K) == 2)
      if (K(1) < 0 | K(2) < 0)
		   mssim = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   mssim = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

if (nargin == 5)
   [H W] = size(window);
   if ((H*W) < 4 | (H > M) | (W > N))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   if (length(K) == 2)
      if (K(1) < 0 | K(2) < 0)
		   mssim = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   mssim = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

C1 = (K(1)*L)^2;
C2 = (K(2)*L)^2;
window = window/sum(sum(window));
img1 = double(img1);
img2 = double(img2);

mu1   = filter2(window, img1, 'valid');
mu2   = filter2(window, img2, 'valid');
mu1_sq = mu1.*mu1;
mu2_sq = mu2.*mu2;
mu1_mu2 = mu1.*mu2;
sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq;
sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq;
sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2;

if (C1 > 0 & C2 > 0)
   ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2));
else
   numerator1 = 2*mu1_mu2 + C1;
   numerator2 = 2*sigma12 + C2;
	denominator1 = mu1_sq + mu2_sq + C1;
   denominator2 = sigma1_sq + sigma2_sq + C2;
   ssim_map = ones(size(mu1));
   index = (denominator1.*denominator2 > 0);
   ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index));
   index = (denominator1 ~= 0) & (denominator2 == 0);
   ssim_map(index) = numerator1(index)./denominator1(index);
end

mssim = mean2(ssim_map);

return


text:
function x=v_randvec(n,m,c,w,mode)
%V_RANDVEC  Generate real or complex GMM/lognormal random vectors X=(N,M,C,W,MODE)
% generates a random matrix of size (|n|,p) where p is the maximum
% dimension of M or C (see note below about row versus column vectors)
%  Inputs:  N        is the number of points to generate
%           M(K,P)   is the mean vectors (one row per mixture)
%           C(K,P)   are diagonal covariances (one row per mixture)
%        or C(P,P,K) are full covariance matrices (one per mixture)
%           W(K)     are the mixture weights (or omit if all mixtures have equal weight)
%           MODE     character string specifying options:
%                       g = real gaussian [default]
%                       c = complex gaussian
%                       l = lognormal
%
% Outputs:  X(N,P) is the output data
%
% Note this routine generates row vectors such that E((x-m)'*(x-m)) = C = cov(x). If
% Alternatively x=(n,m',c,w,mode)' will generate column vectors satisfying E((x-m)*(x-m)') = C = cov(x').

% Bugs/suggestions
%    (1)  New mode 'x' to approximate chi squared

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: v_randvec.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% first sort out the input arguments

sm=size(m);
if nargin<3
    c=ones(sm);         % default to unit variance
end
sc=size(c);
p=max(sm(2),sc(2));     % data dimension
k=sm(1);                % number of mixtures
fullc=(length(sc)>2) || (sc(1)>k);
if nargin<4
    mode='g';   % default to gaussian
    w=ones(k,1);
else
    if ischar(w)
        mode = w;       % w argument has been omitted
        w=ones(k,1);
    elseif nargin<5
        mode='g';
    end
end
ty=mode(1);   % ignore all but first character for type
x=zeros(n,p);   % initialize output array
if sm(2)~=p
    m=repmat(m,1,p);    % if m is given as a scalar
end
if sc(2) ~=p
    c=repmat(c,1,p);    % if c is given as a scalar
end
q=sqrt(0.5);
if k>1
    kx=v_randiscr(w,n);
else
    kx=ones(n,1);
end
for kk=1:k
    nx=find(kx==kk);
    nn=length(nx);
    if nn       % check if we need to generate any from mixture kk

        % extract the mean and cov for this mixture

        mm=m(kk,:);     % mean vector
        if fullc        % full covariance matrix
            cc=c(:,:,kk);
            if ty=='l'      % lognormal distribution - convert mean and covariance
                cc=log(1+cc./(mm.'*mm));
                mm=log(mm)-0.5*diag(cc).';
            end
        else
            cc=c(kk,:);
            if ty=='l'      % lognormal distribution - convert mean and covariance
                cc=log(1+cc(:).'./mm.^2);
                mm=log(mm)-0.5*cc;
            end
        end

        % now generate nn complex or real values

        if ty=='c'  % generate complex values
            xx=q*randn(nn,p)+1i*q*randn(nn,p); % complex-valued unit variance values
        else
            xx=randn(nn,p);   % real-valued unit variance values
        end;

        % scale by the square root of the covariance matrix

        if fullc   % full covariance covariance
            [v,d]=eig((cc+cc')/2);   % force covariance matrix to be hermitian
            xx=(xx.*repmat(sqrt(abs(diag(d))).',nn,1))*v'+repmat(mm,nn,1); % and also positive definite
        else
            xx=xx.*repmat(sqrt(abs(cc)),nn,1)+repmat(mm,nn,1); % different mean/cov for each column
        end
        x(nx,:)=xx;
    end
end
if ty=='l'  % lognormal distribution
    x=exp(x);
end
if ~nargout
    if p==1
        if ty=='c'
            plot(real(x), imag(x),'+');
            xlabel('Real');
            ylabel('Imag');
        else
            nbin=max(min(floor(sqrt(n)),50),5);
            hist(x,nbin);
            xlabel('X');
            ylabel('Frequency');
        end
    else
        [vv,iv]=sort(var(x,0,1));
        iv=sort(iv(end-1:end));
        plot(real(x(:,iv(1))), real(x(:,iv(2))),'+');
        if ty=='c'
            xylab='Real[ x(%d) ]';
        else
            xylab='x(%d)';
        end
        xlabel(sprintf(xylab,iv(1)));
        ylabel(sprintf(xylab,iv(2)));
    end
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [params] = sr_calcSRC(params)
% function [params] = sr_calcZ(params)
% -------------------------------------
% Calculation of seismic rate changes. This function calculates the rates
% for all grid points together.
%
% Input parameters:
% please see function ~/zmap/src/thomas/seismicrates/get_parameter.m
%
% Output parameters:
%   Same as input parameters including but additionally
%   params.mResult_    z,probability of z, beta, probability of beta, and
%                      resolution. This data will be stored into
%                      params.mResult{n}, where n is the number of runs
%                      (see sr_startZ.m)
%   params.mVar        List of all combinations of parameters used for
%                      searching the parameter space ([N, Tw, Tbin]).
%   params.m1          Max. value for z ([z,Resolution,N,Tw,Tbin])
%   params.m2          Max. value for p(z) ([z,Resolution,N,Tw,Tbin])
%   params.m3          Max. value for beta ([z,Resolution,N,Tw,Tbin])
%   params.m4          Max. value for p(beta) ([z,Resolution,N,Tw,Tbin])
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Th. van Stiphout vanstiphout@sed.ethz.ch
% last update: 16.08.2005
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

report_this_filefun(mfilename('fullpath'));

% Initialize
bChk=false;    % for debugging the code

% Determine time period of catalog
params.fTminCat = min(params.mCatalog(:,3));
params.fTmaxCat = max(params.mCatalog(:,3));
% Adjust to decimal years
% fTimePeriod =params.fTimePeriod/365;
% mCat=params.mCatalog(:,3);
mLoc=[params.mCatalog(:,1) params.mCatalog(:,2) params.mCatalog(:,7) ];
% Init result matrix
% mValueGrid_ = [];
%  Selection criteria for subcatalog (between StartTime and
% TimeCut + WindowLength

% Create Indices of the catalog for each Node and select quakes in time
[params.caNodeIndices, params.vResolution] = ex_CreateIndexCatalog(params.mCatalog,...
    params.mPolygon, params.bMap, 0, size(params.mCatalog,1),...
    params.fRadius, params.fSizeRectHorizontal, params.fSizeRectDepth);

% prepare input parameter matrix  [ vN   vTw   vTbin ]
mVar(:,1)=repmat(params.vN,size(params.vTw,1)*size(params.vTbin,1),1);
mVar(:,2)=reshape(repmat(params.vTw',size(params.vN,1)*size(params.vTbin,1),1),...
    size(params.vTw,1)*size(params.vN,1)*size(params.vTbin,1),1);
mVar(:,3)=repmat(reshape(repmat(params.vTbin',size(params.vN,1),1),...
    size(params.vTbin,1)*size(params.vN,1),1),size(params.vTw,1),1);

% loop over all input parameters
for i=1:size(mVar,1)

    % get input parameters
    fTstart=params.fTstart;
    fT=params.fT;
    fTw=mVar(i,2);
    nN=mVar(i,1);
    nTbin=mVar(i,3);
    mNumDeclus_=params.mNumDeclus(:,end);
%     sum(params.mNumDeclus  )
    % vSelMag=(params.mCatalog(:,6)>=params.fMc)

    % create vectors for time periods
    vSelR0_=(params.mCatalog(:,3) > fTstart);
    vSelR1_=(params.mCatalog(:,3) <= fT);
    vSelR2_=(params.mCatalog(:,3) >= fT-fTw);
    vSelR3_=(params.mCatalog(:,3)  <= fT-fTw);

    % combine time periods with decluster info's
    mSelD0_=repmat(vSelR0_,1,size(mNumDeclus_,2)) & mNumDeclus_;
    mSelD1_=repmat(vSelR1_,1,size(mNumDeclus_,2)) & mNumDeclus_;
    mSelD2_=repmat(vSelR2_,1,size(mNumDeclus_,2)) & mNumDeclus_;
    mSelD3_=repmat(vSelR3_,1,size(mNumDeclus_,2)) & mNumDeclus_;

    % only events during certain period
    mSelD00=( mSelD0_ & mSelD1_);   % whole period
    mSelD10=( mSelD0_ & mSelD3_);   % 1st period
    mSelD20=( mSelD1_ & mSelD2_);   % 2nd period

    % events in whole period  (us it for calculating probabilities of z and beta)
    mCat=params.mCatalog(mSelD00,3);

    % repeat matrix for multiplication with grid node distance-ranking
    aSelD02=repmat(mSelD00,1,size(params.mPolygon,1));
    % if gridding is based on radius
    if params.nGriddingMode==1
        for pp=1:size(params.mPolygon,1)
            vNaN(:,pp)=params.vResolution{pp}>params.fRadius;
%             params.caNodeIndices{pp}(vNaN(:,pp))=NaN;
            vResolutionN_(pp)=sum(not(vNaN(:,pp)));
        end
        nN=max(vResolutionN_);
    end
    % transform cells of node indices to matrix (column wise)
    caNodeIndices_=cell2mat(params.caNodeIndices');   % distance-ranking of each event from each grid node

    % create help matrices
    tmp=zeros(size(aSelD02));
    tmp1=(0:1:size(aSelD02,2)-1)*size(params.mCatalog,1);
    tmp2=repmat(tmp1,size(aSelD02,1),1);  %

    % sort rankings for each grid node
    [X1, I1] =sort(caNodeIndices_);
    % do reverse "sort ranking" for each grid node
    [X2, I2] =sort(I1);

    I1=I1+tmp2;I2=I2+tmp2;
    bSelD02=aSelD02(I2); % resort 1-0 matrix
    cSelD02=cumsum(bSelD02);
    dSelD02=logical(bSelD02.*(cSelD02<=nN)); % only the first nN events
    mSelRes=dSelD02;  % copy for the first nN events
    dSelD02=dSelD02(I1); % transform back  / nN Events over whole catalog per grid node
    % transformation back of NaN's-list for event that are not in a certain params.fRadius
    if params.nGriddingMode==1
        vNaN1=vNaN(I1);
    end
    %     eSelD02=bSelD02(I1); % transform back
    %     Checkpoint
    if bChk
        ii=ceil(rand(1)*length(params.mPolygon));
        figure;plot(mLoc(:,1),mLoc(:,2),'.c');
        hold on;plot(params.mPolygon(ii,1),params.mPolygon(ii,2),'ok');
        %  hold on;plot(mLoc(I2(1:nN,ii)-(ii-1)*size(mLoc,1),1),mLoc(I2(1:nN,ii)-(ii-1)*size(mLoc,1),2),'dr')
    end

    % reshape resolution matrix
    vResolution_=cell2mat(params.vResolution');

    % prepare catalog (origin times) - copy them column-wise for each
    % grid point
    mCat0=params.mCatalog(:,3);
    mCat0=(repmat(mCat0,1,size(dSelD02,2)));

    if bChk
        % and for latitude and longitude aswell
        mLon0=params.mCatalog(:,1);
        mLon0=(repmat(mLon0,1,size(dSelD02,2)));
        mLat0=params.mCatalog(:,2);
        mLat0=(repmat(mLat0,1,size(dSelD02,2)));
    end

    % for gridding mode 1 (radius)
    if params.nGriddingMode==1
        adder=repmat(0:size(caNodeIndices_,1):size(caNodeIndices_,1)*(size(caNodeIndices_,2)-1),size(caNodeIndices_,1),1);
        caNodeIndices2_=caNodeIndices_+adder;
        caNodeIndices3_=caNodeIndices2_(~isnan(caNodeIndices2_));
        caNodeIndices4_=zeros(size(caNodeIndices_));
        caNodeIndices4_(caNodeIndices3_)=1;
        caNodeIndices4_=logical(caNodeIndices4_);
%         hold on;plot(params.mCatalog(caNodeIndices4_(:,ii),1),...
%                                  params.mCatalog(caNodeIndices4_(:,ii),2),'sr')
         mCat0(not(caNodeIndices4_))=NaN;
         mLon0(not(caNodeIndices4_))=NaN;
         mLat0(not(caNodeIndices4_))=NaN;

         clear caNodeIndices2_ caNodeIndices3_
    end

    % prepare matrices for subset of catalog
    gSelD02=ones(size(dSelD02))*NaN;
    mCat00=gSelD02; % whole catalog
    mCat10=gSelD02; % 1st period
    mCat20=gSelD02; % 2nd period
    mLon00=gSelD02; mLon00_=mLon00;
    mLat00=gSelD02;    mLat00_=mLat00;
    % put the times in the NaN's matrix
    gSelD02(dSelD02)=(mCat0(dSelD02));

    % whole catalog
    mSelD00=repmat(mSelD00,1,size(dSelD02,2));
    mSelD01=(dSelD02&mSelD00);
    mCat00(mSelD01)=(mCat0(mSelD01));
    if params.nGriddingMode==1
        mCat00(vNaN1)=NaN;
    end
    % Checkpoint
    if bChk
        % figure;histogram(mCat00(:,1),100)
        % Checkpoint plot location
        mLon00=mLon00_;mLat00=mLat00_;
        mLon00(mSelD01)=(mLon0(mSelD01));
        mLat00(mSelD01)=(mLat0(mSelD01));
%         hold on;plot(mLon00(:,ii),mLat00(:,ii),'xk','MarkerSize',10);
        hold on;plot(mLon00(not(vNaN1(:,ii)),ii),mLat00(not(vNaN1(:,ii)),ii),'xk','MarkerSize',10);
    end

    % 1st period
    mSelD10=repmat(mSelD10,1,size(dSelD02,2));
    mSelD11=(dSelD02&mSelD10);
    mCat10(mSelD11)=(mCat0(mSelD11));
    if params.nGriddingMode==1
        mCat10(vNaN1)=NaN;
    end
    % Checkpoint
    if bChk
        % figure;histogram(mCat10(:,1),100)
        % Checkpoint plot location
        mLon10=mLon00_;mLat10=mLat00_;
        mLon10(mSelD11)=(mLon0(mSelD11));
        mLat10(mSelD11)=(mLat0(mSelD11));
%         hold on;plot(mLon10(:,ii),mLat10(:,ii),'dr','MarkerSize',10);
        hold on;plot(mLon10(not(vNaN1(:,ii)),ii),mLat10(not(vNaN1(:,ii)),ii),'dr','MarkerSize',10);
    end

    % 2nd period
    mSelD20=repmat(mSelD20,1,size(dSelD02,2));
    mSelD21=(dSelD02&mSelD20);
    mCat20(mSelD21)=(mCat0(mSelD21));
    if params.nGriddingMode==1
        mCat20(vNaN1)=NaN;
    end
    % Checkpoint
    if bChk
        % figure;histogram(mCat20(:,1),100)
        % Checkpoint plot location
        mLon20=mLon00_;mLat20=mLat00_;
        mLon20(mSelD21)=(mLon0(mSelD21));
        mLat20(mSelD21)=(mLat0(mSelD21));
%         hold on;plot(mLon20(:,ii),mLat20(:,ii),'sb','MarkerSize',10);
        hold on;plot(mLon20(not(vNaN1(:,ii)),ii),mLat20(not(vNaN1(:,ii)),ii),'sb','MarkerSize',10);
    end

    %  lta
    [mLTA, mLTAprob] =calc_zlta(mCat,mCat00,mCat20,params.fTstart, fT,...
        fTw,nTbin, nN);
    % mBeta
    [mBeta, mBetaprob] =calc_beta(mCat,mCat00,mCat20,params.fTstart, fT,...
        fTw,nTbin, nN);
    %      mBeta(:,kk)=calc_beta(mCat10,mCat20,params.fStartTime,...
    %          params.fTimeCut,params.fTwLength,params.fTimeSteps);

    if params.nGriddingMode==1
        % Resolution by No.of Events in a certain distance from grid node,
        % defined by params.fRadius .
        mResolution=vResolutionN_';
    else
        % Resolution by N-th event in radii-vector
        mResolution=max(vResolution_.*mSelRes)';
    end
    % preallocation of mResult_
    if i==1
        mResult_=ones(size(params.mPolygon,1),5,params.nMCS).*NaN;
    end
    mResult_(:,:,i)=[mLTA mLTAprob mBeta mBetaprob double(mResolution)];

    sString=sprintf('End loop %d with N=  %d  Tw=  %d   Tbin=  %d   ',i,...
        mVar(i,1),mVar(i,2),mVar(i,3));
    disp(sString);
end % end loop variation over input parameter

% params.mResult_=mResult_;
% saving parameter space
% params.mVar=mVar;

% Maxima search
[mMaxValue, mMaxIndices]=max(mResult_,[],3);
% ask for parameter settings (parameter settings for max for each  grid
% point)
% nMax=2; % (i.e. 1: z, 2:prob(z), 3:beta, 4:prob(beta),5:radius)
mParamax1=mVar(mMaxIndices(:,1),1:3);
mParamax2=mVar(mMaxIndices(:,2),1:3);
mParamax3=mVar(mMaxIndices(:,3),1:3);
mParamax4=mVar(mMaxIndices(:,4),1:3);

% plotting preliminary results / data
if bChk
    % plot results for grid points with contour lines that indicate z,
    % prob(z), beta, prob(beta)
    %     figure;pcolor( params.mPolygon(:,1), params.mPolygon(:,2), mResult_(:,1))
    figure;
    XX_=repmat(params.vN,1,size(params.vTw,1));
    YY_=repmat(params.vTw,1,size(params.vN,1))';
    for ii_=1:10
        ngp_=ceil(rand(1,1)*1170);
        for kk_=1:4
            subplot(2,2,kk_);
            contour(XX_,YY_,reshape(squeeze(mResult_(ngp_,kk_,:)),...
                size(params.vN,1),size(params.vTw,1)));
            colorbar;
        end
        pause
    end
    %     pcolor(reshape(squeeze(mResult_(1,1,:)),size(params.vN,1),size(params.vTw,1)));
end

% Prepare mResult_
% vSel=find(mResult_(:,5) > params.fMaxRadius);
% mResult_(vSel,1:4)=NaN';

% Pre Result matrix (result for 1
params.m1=[mMaxValue(:,1) double(mResolution) mParamax1];
params.m2=[mMaxValue(:,2) double(mResolution) mParamax2];
params.m3=[mMaxValue(:,3) double(mResolution) mParamax3];
params.m4=[mMaxValue(:,4) double(mResolution) mParamax4];
% figure;pcolor(reshape(params.mPreResult(:,1),39,30));
% colorbar

author: CelsoReyes


text:
function centerGUI(h)

OldUnits = get(h,'Units');
set(h,'Units','pixels');
OldPos = get(h,'Position');
FigWidth = OldPos(3);
FigHeight = OldPos(4);

ScreenUnits = get(0,'Units');
set(0,'Units','pixels');
ScreenSize = get(0,'MonitorPosition'); % MonitorPosition works better than ScreenSize for multiheaded comps
ScreenSize = ScreenSize(1,:);
set(0,'Units',ScreenUnits);

FigPos(1)=1/2*(ScreenSize(3)-FigWidth);
FigPos(2)=2/3*(ScreenSize(4)-FigHeight);
FigPos(3:4)=[FigWidth FigHeight];
set(h,'Position',FigPos);
set(h,'Units',OldUnits);

end


text:
function varargout=randfilt(varargin)
%
% For calling details please see v_randfilt.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_randfilt.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: randfilt.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_randfilt(varargin{:});
else
    v_randfilt(varargin{:});
end


text:
function [fragments, junctions, neighbor_data, seg, avg_seg_colors, ...
    polyfragments, poly_params] = seg2fragments(seg, img, min_area, order)
%
% [fragments, junctions, neighbor_data, new_seg, avg_seg_colors, ...
%  <polyfragments>, <poly_params>] = seg2fragments(seg, img, min_area, <order>)
%
%  Finds the boundary fragments for an oversegmentation of an image:
%
%  (1) Takes in an over-segmentation labeling (with labels 1:N),
%  (2) Fills any pixels with a label of zero with the neighboring label
%      that has the most similar color value in the original RGB color
%      image 'img',
%  (3) Cleans up segments which have a single pixel jutting diagonally off
%      of a corner somewhere, by reassigning these offending pixels to the
%      neighboring segment with the most similar color (as in (2)).
%  (4) Breaks diagonal-only connected segments into two new segments.
%  (5) Merges segments that are too small with a neighbor with the most
%      similar average color, and
%  (6) Chains fragments along the boundaries between segments, stopping at
%      junctions and image borders.
%
%  Returns the fragments as a cell array of [Mx2] matrices which contain
%  the M (x,y) coordinates of of each fragment's constituent elements, the
%  fragment neighbor information (see below for details), and the new
%  (filled-in / merged) segmentation.
%
%  'neighbor_data' is a struct with two fields:
%   - 'junction_fragmentlist': a list which contains for each junction an
%                              array of indices of the fragments that meet
%                              there
%   - 'fragment_junctionlist': a list which contains for each fragment the
%                              indices of the junctions which bound it
%   - 'fragment_segments': a list which contains for each fragment the ID
%                          numbers of the segments found to the left and
%                          right of that fragment, stored [left right].
%   - 'segment_fragments': a list which contains for each segment the index
%                          of the fragments that border it
%
%  Thus, the indices of fragment i's neighbors on one end can be found by:
%    junction_fragmentlist{fragment_junctionlist{i}(1)}
%
%  If requested, also fits polynomials of degree 'order' to each fragment
%  and returns those as well (and their parameters).
%

% DEPENDENCIES:
%  - RGB2Lab.m
%  - cracks2fragments.m
%  - fill_in_segmentation.m
%  - fit_curves.m
%    - fit_poly_to_fragment.m
%  - remove_small_segments.m
%  - seg2cracks.m
% [- drawedges.m (will use it if available, not a problem if not)]

[nrows, ncols] = size(seg);

% Convert the image to Lab space so that color distances used below are
% more meaningful
if ~isempty(img)
    img_lab = RGB2Lab(img);
end

while (any(seg(:)==0))
    seg = fill_in_segmentation(img, seg, 0, 4);
end

num_segments = max(seg(:));

% Find locations where there's a pixel jutting off diagonally from  a
% corner of a segment.  We want to remove those and attach them to the
% nearest-colored 4-connected neighbor.
% Note: this is actually redundant since this would also be covered below,
% but it's much faster to remove the simple singleton stragglers like this
% first (since it requires fewer calls to bwselect).
% stragglers = find( seg~=seg(:,[2:end end]) & seg~=seg(:,[1 1:end-1]) & ...
%     seg~=seg([2:end end],:) & seg~=seg([1 1:end-1],:) );
% if(any(stragglers(:)))
%     corner_neighbor_offsets = [-1-nrows -1+nrows 1-nrows 1+nrows];
%
%     for(i=1:length(stragglers))
%         % if the straggler was diagonnally connecting >=2 larger pieces of
%         % this segment, we need to re-label all but one of them to avoid
%         % inadvertantly splitting a segment into multiple pieces by
%         % removing the straggler
%         [y,x] = ind2sub([nrows ncols], stragglers(i));
%         corner_x = x + [-1 1 1 -1];
%         corner_y = y + [-1 -1 1 1];
%         out_of_bounds = find(corner_x<1 | corner_x>ncols | corner_y<1 | corner_y>nrows);
%         corner_x(out_of_bounds) = [];
%         corner_y(out_of_bounds) = [];
%         corner_neighbors = sub2ind([nrows ncols], corner_y, corner_x);
%         to_relabel = find( seg(corner_neighbors)==seg(stragglers(i)) );
%         if(length(to_relabel)>=2)
%             bw_img = seg==seg(stragglers(i));
%             for(j = 1:length(to_relabel))
%                 num_segments = num_segments+1;
%                 seg(bwselect(bw_img,corner_x(to_relabel(j)),corner_y(to_relabel(j)),4)) = num_segments;
%             end
%         end
%     end
%
%     % Now actually get rid of the stragglers by setting them to zero and
%     % filling them in:
%     seg(stragglers) = 0;
%     seg = fill_in_segmentation(img, seg, 0, 4);
%
% end

% Convert the segmentation to a crack image.  Note that we need this for
% the next step, even if we have to end up re-doing it later if small
% segments get removed.
cracks = seg2cracks(seg);

% Find locations where a single segment narrows to a singe-pixel diagonal
% connection.  We want to split such locations into two segments.
down_right_splits = find(cracks==15 & seg==seg([2:end end],[2:end end]));
down_left_splits  = find(cracks==15 & seg(:,[2:end end])==seg([2:end end],:));
[r,c] = ind2sub([nrows ncols], down_right_splits);
for(i=1:length(down_right_splits))
    % Need to set _both_ of the resulting segments after the split to new
    % ID numbers to avoid inadvertantly forming a non-contiguous segment
    % later... (argh, twice as many bwselect calls!)
    temp = seg==seg(down_right_splits(i));
    num_segments = num_segments+1;
    seg(bwselect(temp,c(i),r(i),4)) = num_segments;
    num_segments = num_segments+1;
    seg(bwselect(temp,c(i)+1,r(i)+1,4)) = num_segments;
end
[r,c] = ind2sub([nrows ncols], down_left_splits);
for(i=1:length(down_left_splits))
    temp = seg==seg(r(i),c(i)+1);
    num_segments = num_segments+1;
    seg(bwselect(temp,c(i)+1,r(i),4)) = num_segments;
    num_segments = num_segments+1;
    seg(bwselect(temp,c(i),r(i)+1,4)) = num_segments;
end

if(nargin < 3 || isempty(min_area))
    min_area = 0;
end

if min_area > 0
    [seg, avg_seg_colors] = remove_small_segments(seg, img_lab, min_area);
end

% Assuming any segments were removed, we'll need to update the crack info
cracks = seg2cracks(seg);

try
    [fragments, junctions, neighbor_data] = cracks2fragments(cracks, seg);
catch
    disp('warning: no isolation check');
    [fragments, junctions, neighbor_data] = cracks2fragments(cracks, seg, 0);
end

if(nargout>=5 || nargout==0)
    if(nargin < 4)
        order = 3;
    end

    [polyfragments, poly_params] = fit_curves(fragments, order);
end

if(nargout==0)
    figure(gcf)
    subplot 121, hold off
    imagesc(img), axis image, hold on
    title('Segment Borders')

    subplot 122, hold off
    imagesc(img), axis image, hold on
    title(['Polynomial Fits to Fragments (order=' num2str(order)])

    if(exist('file', 'drawedges'))
        subplot 121, drawedges(fragments, 'rand');
        subplot 122, drawedges(polyfragments, 'rand');
    else
        subplot 121
        for(i=1:length(fragments))
            plot(fragments{i}(:,1), fragments{i}(:,2), 'r', 'LineWidth', 2);
        end
        subplot 122
        for(i=1:length(polyfragments))
            plot(polyfragments{i}(:,1), polyfragments{i}(:,2), 'r', 'LineWidth', 2);
        end
    end
end

author: Cloud-CV

repo: object-proposals


text:
%Example of using MUMPS in matlab with sparse right hansd side

% initialization of a matlab MUMPS structure
id = initmumps;
id = dmumps(id);
load lhr01;
mat = Problem.A;
% JOB = 6 means analysis+facto+solve
id.JOB = 6;
% we set the rigth hand side
id.RHS = ones(size(mat,1),2);
id.RHS(:,2) = 2*id.RHS(:,2);
id.RHS = sparse(id.RHS);
%call to mumps
id = dmumps(id,mat);
if(norm(mat*id.SOL - id.RHS,'inf') > sqrt(eps))
	disp('WARNING : precision may not be OK');
else
	disp('SOLUTION OK');
end
norm(mat*id.SOL - id.RHS,'inf')
% solution OK
% destroy mumps instance
id.JOB = -2;
id = dmumps(id)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
classdef dme_gen3p_opf < mp.dme_gen3p & mp.dme_shared_opf

%   MATPOWER
%   Copyright (c) 2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end     %% properties

%     methods
%     end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function plotQualityMeasuresForCompressedDatasets(srResults, qualityMeasure, binningFactors, compression, srMethods, normalized, outputDir)

    figure_bin2 = plotQualityMeasuresAsBarplot(srResults, qualityMeasure, binningFactors(1), srMethods, compression, normalized);
    figure_bin3 = plotQualityMeasuresAsBarplot(srResults, qualityMeasure, binningFactors(2), srMethods, compression, normalized);
    figure_bin4 = plotQualityMeasuresAsBarplot(srResults, qualityMeasure, binningFactors(3), srMethods, compression, normalized);
    LinkFigures([figure_bin2 figure_bin3 figure_bin4]);

    if normalized
        qualityMeasure = [qualityMeasure, '_norm'];
    end

    % Save to Tikz.
    figure(figure_bin2);
    matlab2tikz([outputDir, '/', 'compressionDatasets_bin2_', qualityMeasure, '_bars.tikz'], ...
        'height', '\figureheight', 'width', '\figurewidth', 'showInfo', false, ...
        'extraaxisoptions', ['xlabel near ticks,', ...
        'ylabel near ticks,', ...
        'scaled y ticks=false,', ...
        'yticklabel style={/pgf/number format/fixed, /pgf/number format/precision=2},']);

    figure(figure_bin3);
    matlab2tikz([outputDir, '/', 'compressionDatasets_bin3_', qualityMeasure, '_bars.tikz'], ...
        'height', '\figureheight', 'width', '\figurewidth', 'showInfo', false, ...
        'extraaxisoptions', ['xlabel near ticks,', ...
        'ylabel near ticks,', ...
        'scaled y ticks=false,', ...
        'yticklabel style={/pgf/number format/fixed, /pgf/number format/precision=2},']);

    figure(figure_bin4);
    matlab2tikz([outputDir, '/', 'compressionDatasets_bin4_', qualityMeasure, '_bars.tikz'], ...
        'height', '\figureheight', 'width', '\figurewidth', 'showInfo', false, ...
        'extraaxisoptions', ['xlabel near ticks,', ...
        'ylabel near ticks,', ...
        'scaled y ticks=false,', ...
        'yticklabel style={/pgf/number format/fixed, /pgf/number format/precision=2},']);

function figureNumber = plotQualityMeasuresAsBarplot(srResults, qualityMeasure, binningFactor, srMethods, compression, normalized)

    % Compute mean of the quality measure for different SR methods and
    % binning factors.
    barMeans = zeros(length(compression), length(srMethods));
    for compressionLevelIdx = 1:length(compression)
        for srMethodIdx = 1:length(srMethods)
            % Get quality measures for given SR method and binning factor.
            qualityMeasureVals = [srResults([ strcmp({srResults.compressionLevel}, compression{compressionLevelIdx}) ] ...
                & [srResults(:).binningFactor] == binningFactor ...
                &  strcmp(srMethods{srMethodIdx}, {srResults(:).srMethod}) > 0).(qualityMeasure)];

            if normalized
                % Normalize the quality measure with the input
                % low-resolution image (after NN interpolation).
                qualityMeasureVals_nn = [srResults([ strcmp({srResults.compressionLevel}, compression{compressionLevelIdx}) ] ...
                    & [srResults(:).binningFactor] == binningFactor &  strcmp('nn', {srResults(:).srMethod}) > 0).(qualityMeasure)];
                qualityMeasureVals = (qualityMeasureVals - qualityMeasureVals_nn) ./ abs(qualityMeasureVals_nn);
            end

            % Compute the mean of the quality measure.
            barMeans(compressionLevelIdx, srMethodIdx) = mean(qualityMeasureVals);
        end
    end

    % Plot mean quality measures as bar plot.
    figure;
    figureNumber = get(gcf, 'Number');
    bar_handle = bar(barMeans);
    for k = 1:length(srMethods)
        set(bar_handle(k), 'FaceColor', getFaceColorForSRMethod(srMethods{k}));
    end
    hold off;
    pbaspect([4 1.25 1]);
    set(gca, 'XTickLabel', {'Uncoded', 'QP10', 'QP20', 'QP30', 'QP40'}, 'XTick', 1:(5*length(srMethods)));
    ymax = max(1.05*max(barMeans(:)));
    ymin = max(0.99*min(barMeans(:)), -0.25*ymax);
    ylim([ymin ymax]);
    ylabelStr = upper(qualityMeasure);
    if normalized
        ylabelStr = ['Norm. ', ylabelStr]; %#ok<AGROW>
    end
    ylabel(ylabelStr);
    grid on;

author: thomas-koehler

repo: SupER


text:
function [opts, args] = vl_argparse(opts, args, varargin)
%VL_ARGPARSE Parse list of parameter-value pairs.
%   OPTS = VL_ARGPARSE(OPTS, ARGS) updates the structure OPTS based on
%   the specified parameter-value pairs ARGS={PAR1, VAL1, ... PARN,
%   VALN}. If a parameter PAR cannot be matched to any of the fields
%   in OPTS, the function generates an error.
%
%   Parameters that have a struct value in OPTS are processed
%   recursively, updating the individual subfields.  This behaviour
%   can be suppressed by using VL_ARGPARSE(OPTS, ARGS, 'nonrecursive'),
%   in which case the struct value is copied directly (hence deleting any
%   existing subfield existing in OPTS). A direct copy occurrs also if the
%   struct value in OPTS is a structure with no fields. The nonrecursive
%   mode is especially useful when a processing time is a concern.
%
%   One or more of the (PAR, VAL) pairs in the argument list can be
%   replaced by a structure; in this case, the fields of the structure
%   are used as paramater names and the field values as parameter
%   values. This behaviour, while orthogonal to structure-valued parameters,
%   is also disabled in the 'nonrecursive' mode.
%
%   A shorthand notation for substructs is available: the pair
%   ('field.subfield', value) creates a substruct, i.e. ('field',
%   struct('subfield', value)). It also works recursively (sub-sub-fields).
%
%   [OPTS, ARGS] = VL_ARGPARSE(OPTS, ARGS) copies any parameter in
%   ARGS that does not match OPTS back to ARGS instead of producing an
%   error. Options specified as structures are passed back as a list
%   of (PAR, VAL) pairs.
%
%   A further option is VL_ARGPARSE(OPTS, ARGS, 'merge'), which merges any
%   parameter in ARGS that do not match OPTS into OPTS. This essentially
%   turns off errors for unknown options, merging them anyway, and is
%   useful if OPTS does not necessarily define all possible default values.
%   The second output argument will always be empty. Another view is that
%   VL_ARGPARSE(OPTS1, OPTS2, 'merge') merges two structs into one.
%
%   Example::
%     The function can be used to parse a list of arguments
%     passed to a MATLAB functions:
%
%        function myFunction(x,y,z,varargin)
%        opts.parameterName = defaultValue ;
%        opts = vl_argparse(opts, varargin)
%
%     If only a subset of the options should be parsed, for example
%     because the other options are interpreted by a subroutine, then
%     use the form
%
%        [opts, varargin] = vl_argparse(opts, varargin)
%
%     that copies back to VARARGIN any unknown parameter.
%
%   See also: VL_HELP().

% Copyright (C) 2015-16 Andrea Vedaldi and Karel Lenc.
% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% Tishis file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if ~isstruct(opts) && ~isobject(opts), error('OPTS must be a structure') ; end
if ~iscell(args), args = {args} ; end

recursive = true ;
merge = false ;

if numel(varargin) > 2
  error('There can be at most two options.') ;
end

for i = 1:numel(varargin)
  switch lower(varargin{i})
    case 'nonrecursive'
      recursive = false ;
    case 'merge'
      merge = true ;
    otherwise
      error('Unknown option specified.') ;
  end
end

optNames = fieldnames(opts)' ;

% convert ARGS into a structure
ai = 1 ;
keep = false(size(args)) ;
while ai <= numel(args)

  % Check whether the argument is a (param,value) pair or a structure.
  if recursive && isstruct(args{ai})
    params = fieldnames(args{ai})' ;
    values = struct2cell(args{ai})' ;
    if nargout <= 1
      opts = vl_argparse(opts, vertcat(params,values), varargin{:}) ;
    else
      [opts, rest] = vl_argparse(opts, reshape(vertcat(params,values), 1, []), varargin{:}) ;
      args{ai} = cell2struct(rest(2:2:end), rest(1:2:end), 2) ;
      keep(ai) = true ;
    end
    ai = ai + 1 ;
    continue ;
  end

  if ~isstr(args{ai})
    error('Expected either a param-value pair or a structure.') ;
  end

  param = args{ai} ;
  value = args{ai+1} ;

  if any(param == '.')
    % Handy notation for substructs: convert {'train.learningRate', value}
    % to {'train', struct('learningRate', value)}. Nested structs also work
    parts = strsplit(param, '.') ;
    subs = struct('type', repmat({'.'}, 1, numel(parts) - 1), 'subs', parts(2:end)) ;
    param = parts{1} ;
    value = subsasgn(struct(), subs, value) ;
  end

  p = find(strcmpi(param, optNames)) ;
  if numel(p) ~= 1
    if merge  % Merge unknown fields, instead of erroring or leaving it
      field = param ;
    elseif nargout <= 1
      error('Unknown parameter ''%s''', param) ;
    else
      keep([ai,ai+1]) = true ;
      ai = ai + 2 ;
      continue ;
    end
  else
    field = optNames{p} ;
  end

  if ~recursive
    opts.(field) = value ;
  else
    if isfield(opts, field) && isstruct(opts.(field)) && numel(fieldnames(opts.(field))) > 0
      % The parameter has a  non-empty struct value in OPTS:
      % process recursively.
      if ~isstruct(value)
        error('Cannot assign a non-struct value to the struct parameter ''%s''.', ...
          field) ;
      end
      if nargout > 1
        [opts.(field), args{ai+1}] = vl_argparse(opts.(field), value, varargin{:}) ;
      else
        opts.(field) = vl_argparse(opts.(field), value, varargin{:}) ;
      end
    else
      % The parameter does not have a struct value in OPTS: copy as is.
      opts.(field) = value ;
    end
  end

  ai = ai + 2 ;
end

args = args(keep) ;

repo: LapSRN


categnames={'aeroplane';'bicycle';'bird';'boat';'bottle';'bus';'car';'cat';'chair';'cow';'diningtable';'dog';'horse';'motorbike';'person';'pottedplant';'sheep';'sofa';'train';'tvmonitor'};
transport=[1 2 4 6 7 14 19];
artic=[3 8 10 12 13 15 17];
indoor=[5 9 11 16 18 20];

groups={transport, artic, indoor};
for i=1:numel(groups)
for k=1:numel(groups{i})
	similar{groups{i}(k)}=setdiff(groups{i}, groups{i}(k));
end
end
for i=1:20
other{i}=setdiff([1:20], similar{i});
end


text:
function [er] = mountain_car_policy_expected_return(nfq_Q,X)
%MOUNTAIN_CAR_POLICY_EXPECTED_RETURN Computes the average expected return
% of a policy given a Q-value function and set set of Initial starting
% positions , X
%
%   input ------------------------------------
%
%       o Q: either handle or cell array
%
%       o X: (M x D), set of initial starting positions
%

maxsteps = 500;
grafic   = false;
bRecord  = false;
epsilon  = 0;

policy      = @(x,epsilon)nfq_policy(x,nfq_Q,BuildActionList(),epsilon);
reward      = @(x)nfq_mountain_car_reward( x );
start_type  = 'semi-random';

N = size(X,1);

R = zeros(N,1);

for i=1:N
    %                                            maxsteps,epsilon,policy,reward,grafic,bRecord
    disp(['  i(' num2str(i) ')']);
    x_init = X(i,:);
    [~,steps] = car_run_policy(maxsteps,epsilon,policy,reward,grafic,start_type,bRecord,x_init);
    R(i) = steps./maxsteps; % undiscounted total reward
end

er = [mean(R),var(R)];

end

author: epfl-lasa


text:
function aedat = TrimTime(aedat, startTime, endTime, reZero)

%{
default is not to rezero
Times come in as seconds
%}

dbstop if error

if ~exist('startTime', 'var')
    startTime = 0;
end

if ~exist('endTime', 'var') || endTime == 0
    endTime = double(aedat.info.lastTimeStamp) / 1e6;
end

% Choose the right data type to avoid converting the timestamps when a
% comparison is necessary.
if ~isfield(aedat.info, 'fileFormat') || aedat.info.fileFormat == 2
    startTime = uint32(startTime * 1e6);
    endTime = uint32(endTime * 1e6);
else
    startTime = uint64(startTime * 1e6);
    endTime = uint64(endTime * 1e6);
end

if ~isfield(aedat, 'data')
    disp('No data found')
    return
end

%% Special

if isfield(aedat.data, 'special')
    keepLogical = aedat.data.special.timeStamp >= startTime & aedat.data.special.timeStamp <= endTime;
    aedat.data.special.timeStamp = aedat.data.special.timeStamp(keepLogical);
    aedat.data.special.address   = aedat.data.special.address  (keepLogical);
end

%% Polarity

if isfield(aedat.data, 'polarity')
    keepLogical = aedat.data.polarity.timeStamp >= startTime & aedat.data.polarity.timeStamp <= endTime;
    aedat.data.polarity.timeStamp = aedat.data.polarity.timeStamp (keepLogical);
    aedat.data.polarity.x           = aedat.data.polarity.x       (keepLogical);
    aedat.data.polarity.y           = aedat.data.polarity.y       (keepLogical);
    aedat.data.polarity.polarity    = aedat.data.polarity.polarity(keepLogical);
end

%% Frames

if isfield(aedat.data, 'frame')
    keepLogical = aedat.data.frame.timeStampStart >= startTime & aedat.data.frame.timeStampStart <= endTime;
    aedat.data.frame.timeStampStart = aedat.data.frame.timeStampStart(keepLogical);
    aedat.data.frame.timeStampEnd   = aedat.data.frame.timeStampEnd  (keepLogical);
    aedat.data.frame.samples        = aedat.data.frame.samples       (keepLogical);
    aedat.data.frame.xLength        = aedat.data.frame.xLength       (keepLogical);
    aedat.data.frame.yLength        = aedat.data.frame.yLength       (keepLogical);
    aedat.data.frame.xPosition      = aedat.data.frame.xPosition     (keepLogical);
    aedat.data.frame.yPosition      = aedat.data.frame.yPosition     (keepLogical);
end

%% Imu6

if isfield(aedat.data, 'imu6')
    keepLogical = aedat.data.imu6.timeStamp >= startTime & aedat.data.imu6.timeStamp <= endTime;
    aedat.data.imu6.timeStamp = aedat.data.imu6.timeStamp (keepLogical);
    aedat.data.imu6.accelX    = aedat.data.imu6.accelX    (keepLogical);
    aedat.data.imu6.accelY    = aedat.data.imu6.accelY    (keepLogical);
    aedat.data.imu6.accelZ    = aedat.data.imu6.accelZ    (keepLogical);
    aedat.data.imu6.gyroX     = aedat.data.imu6.gyroX     (keepLogical);
    aedat.data.imu6.gyroY     = aedat.data.imu6.gyroY     (keepLogical);
    aedat.data.imu6.gyroZ     = aedat.data.imu6.gyroZ     (keepLogical);
    aedat.data.imu6.temperature = aedat.data.imu6.temperature(keepLogical);
end

%% SAMPLE, EAR TODO

%% Point1D

if isfield(aedat.data, 'point1D')
    keepLogical = aedat.data.point1D.timeStamp >= startTime ...
                & aedat.data.point1D.timeStamp <= endTime;
    aedat.data.point1D.timeStamp = aedat.data.point1D.timeStamp(keepLogical);
    aedat.data.point1D.type      = aedat.data.point1D.type     (keepLogical);
    aedat.data.point1D.x         = aedat.data.point1D.x        (keepLogical);
end

%% Point2D

if isfield(aedat.data, 'point2D')
    keepLogical = aedat.data.point2D.timeStamp >= startTime ...
                & aedat.data.point2D.timeStamp <= endTime;
    aedat.data.point2D.timeStamp = aedat.data.point2D.timeStamp(keepLogical);
    aedat.data.point2D.type      = aedat.data.point2D.type     (keepLogical);
    aedat.data.point2D.x         = aedat.data.point2D.x        (keepLogical);
    aedat.data.point2D.y         = aedat.data.point2D.y        (keepLogical);
end

%% Point3D

if isfield(aedat.data, 'point3D')
    keepLogical = aedat.data.point3D.timeStamp >= startTime ...
                & aedat.data.point3D.timeStamp <= endTime;
    aedat.data.point3D.timeStamp = aedat.data.point3D.timeStamp(keepLogical);
    aedat.data.point3D.type      = aedat.data.point3D.type     (keepLogical);
    aedat.data.point3D.x         = aedat.data.point3D.x        (keepLogical);
    aedat.data.point3D.y         = aedat.data.point3D.y        (keepLogical);
    aedat.data.point3D.z         = aedat.data.point3D.z        (keepLogical);
end

%% Tidy up

aedat = NumEventsByType(aedat);
aedat = FindFirstAndLastTimeStamps(aedat);

%% Rezero

if exist('reZero', 'var') && reZero
    aedat = ZeroTime(aedat);
end

author: panpanfei


text:
function [ishappy, cutoff] = standardCheck(f, values, data, pref)
%STANDARDCHECK   Attempt to trim trailing Fourier coefficients in a TRIGTECH.
%   [ISHAPPY, CUTOFF] = STANDARDCHECK(F) uses the routine STANDARDCHOP to
%   compute a positive integer CUTOFF which represents the number of
%   coefficients of F that are deemed accurate enough to keep.  ISHAPPY is TRUE
%   if the CUTOFF value returned by STANDARDCHOP is less than LENGTH(F) and
%   FALSE otherwise.
%
%   [ISHAPPY, CUTOFF] = STANDARDCHECK(F, VALUES, DATA, PREF) allows additional
%   preferences to be passed. VALUES is a matrix of the function values of F at
%   the corresponding interpolation points. DATA.VSCALE is an approximation of
%   the maximum function value of F on a possibly larger approximation
%   interval.  PREF is a data structure used to pass in additional information,
%   e.g. a target accuracy tolerance could be passed using PREF.CHEBFUNEPS.
%
% See also CLASSICCHECK, STRICTCHECK, LOOSECHECK.

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

% Grab the coefficients of F.
coeffs = abs(f.coeffs(end:-1:1,:));
[n, m] = size(coeffs);

% Compute some values if none were given.
if ( nargin < 2 || isempty(values) )
    values = f.coeffs2vals(f.coeffs);
end

% In order to work with STANDARDCHOP, the coefficients of F are modified so that
% the entries corresponding to wave numbers k and -k appear sequentially in the
% new matrix of coefficients. These entries are also replaced by the sum of the
% absolute values of the k and -k coefficients.

% Need to handle odd/even cases separately.
isEven = mod(n, 2) == 0;
if ( isEven )
    coeffs = [coeffs(n,:) ; coeffs(n-1:-1:n/2+1,:) + coeffs(1:n/2-1,:) ; coeffs(n/2,:)];
else
    coeffs = [coeffs(n:-1:(n+1)/2+1,:) + coeffs(1:(n+1)/2-1,:) ; coeffs((n+1)/2,:)];
end
coeffs = flipud(coeffs);
coeffs = [coeffs(1,:) ; kron(coeffs(2:end,:), [1 ; 1])];

% Initialize ISHAPPY.
ishappy = false;

% Initialize CUTOFF.
cutoff = n;

% NaNs are not allowed.
if ( any(isnan(coeffs)) )
    error('CHEBFUN:TRIGTECH:standardCheck:nanEval', ...
        'Function returned NaN when evaluated.')
end

% Grab some preferences.
if ( nargin == 1 )
    pref = f.techPref();
    tol = pref.chebfuneps;
elseif ( isnumeric(pref) )
    tol = pref;
else
    tol = pref.chebfuneps;
end

% Reshape TOL.
if ( size(tol, 2) ~= m )
  tol = ones(1, m)*max(tol);
end

% Scale TOL by VSCL/||F||;
nrmf = max(abs(values), [], 1);
tol = tol.*data.vscale./nrmf;

% Loop through columns of coeffs
ishappy = false(1, m);
cutoff = zeros(1, m);
for k = 1:m

    % Call STANDARDCHOP.
    cutoff(k) = standardChop(coeffs(:,k), tol(k));

    % Check for happiness.
    ishappy(k) = ( cutoff(k) < n );

    % Divide CUTOFF by 2.
    if ( mod(cutoff(k), 2) == 0 )
        cutoff(k) = cutoff(k)/2;
    else
        cutoff(k) = (cutoff(k) - 1)/2;
    end

    % Break if unhappy.
    if ( ~ishappy(k) )
        break
    end

end

% Set outputs.
ishappy = all(ishappy);

% CUTOFF is always odd.
cutoff = 2*max(cutoff) + 1;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
clear all
close all

Nobs=40;

[M,U] = mci_approach_struct (Nobs);

rand_params=1;
if rand_params
    P = spm_normrnd(M.pE,M.pC,1);
else
    V = 30;
    tau = 8;
    P = [log(V),log(tau)]';
end

yhat = mci_approach_gen (P,M,U);
Y = yhat + sqrt(M.Ce)*randn(Nobs,1);

% Check analytic gradients using finite differences
%dLdp = approach_deriv (P,M,U,Y);
%dLdp = spm_diff(M.L,P,M,U,Y,1);

% figure
% plot(U.X,Y);
% hold on
% plot(U.X,yhat,'r');

%mcmc.inference='vl';
mcmc.inference='langevin';

mcmc.verbose=0;
mcmc.maxits=1024;

post = spm_mci_post (mcmc,M,U,Y,P);

disp('Prior Mean (circle):');
disp(M.pE);

disp('Posterior Mean from MCI (cross)');
disp(post.Ep)

disp('True (plus)');
disp(P)

% Plot posterior surface
S.Nbins=100;
S.pxy(1,:)=linspace(2,4,S.Nbins);
S.pxy(2,:)=linspace(0.5,3.5,S.Nbins);
S.param{1}='P(1)';
S.param{2}='P(2)';
S.name={'log V_a','log \tau'};
%[L,S] = mci_plot_surface (P,M,U,Y,S,'prior');
%[L,S] = mci_plot_surface (P,M,U,Y,S,'like');
[L,S] = mci_plot_surface (P,M,U,Y,S,'post');
hold on
ms=10;
plot(M.pE(1),M.pE(2),'wo','MarkerSize',ms);
plot(post.Ep(1),post.Ep(2),'wx','MarkerSize',ms);
plot(P(1),P(2),'w+','MarkerSize',ms);
j=post.ind;
plot(post.P(1,j),post.P(2,j),'w.');

diag.essplot=1
diag.ind=j;
mess=spm_mci_diag(post,diag);

stats = spm_mci_mvnpost (post,'ESS')
stats = spm_mci_mvnpost (post,'thinning')


text:
function test_suite=BatchExample_qmt_bssfp_test
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 TestSetup
setenv('ISCITEST','1') % go faster! Fit only 2 voxels in FitData.m

function test_batch
curdir = pwd;

tmpDir = tempdir;
mkdir(tmpDir);
cd(tmpDir);

Modellist = {'filter_map'};

for iModel = 1:length(Modellist)
    disp('===============================================================')
    disp(['Testing: ' Modellist{iModel} ' BATCH...'])
    disp('===============================================================')

    eval(['Model = ' Modellist{iModel}]);

    qMRgenBatch(Model,pwd)

    % Test if any dataset exist
    isdata = true;
    try
        Model.onlineData_url;
    catch
        isdata = false;
    end

    % Run Batch
    if isdata
        starttime = tic;
        eval([Modellist{iModel} '_batch'])
        toc(starttime)
    end
    close all
    cd ..

end
cd(curdir)

function TestTeardown
setenv('ISCITEST','0')

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% this m-file calculates the real roots of the given polynomial using
% newton raphson technique.this m-file calls the functions in the two m-files named as syn_division and
% derivate.
% coeff_function is the 1xn matrix conatining the coeff of the polynomial.
% Keerthi Venkateswara Rao
function [final_roots,functionvalue] = newton(coeff_function,initial_guess,error_tolerance,max_iterations)
iterations=0;
max_no_roots=size(coeff_function,2);
final_roots=0;
functionvalue=0;
for no_roots=1:max_no_roots-1
    fun_root_new=initial_guess;
    flag=1;
    coeff_der_function=derivate(coeff_function);
    order_fun=size(coeff_function,2);
    order_der_fun=size(coeff_der_function,2);
    while flag==1
        fun_root_old=fun_root_new;
        fx=0;
        dfx=0;
        nonzero=1;
        while nonzero==1
            powers=order_fun-1;
            for index=1:order_fun
                fx=fx+coeff_function(index)*fun_root_old^powers;
                powers=powers-1;
            end
            powers=order_der_fun-1;
            for index=1:order_der_fun
                dfx=dfx+coeff_der_function(index)*fun_root_old^powers;
                powers=powers-1;
            end
            if dfx==0
                fun_root_old=fun_root_old+1;
            else
                nonzero=0;
            end
        end
        iterations = iterations + 1;
        fun_root_new = fun_root_old - fx/dfx;
        if iterations >= max_iterations
            flag=0;
        elseif  abs(fun_root_new-fun_root_old)<=error_tolerance
            flag=0;
            final_roots(no_roots)=fun_root_new;
            functionvalue(no_roots)=fx;
        end
    end
    coeff_function=syn_division(coeff_function,fun_root_new);
 end


text:
function [A,b,x] = phillips(n)
%PHILLIPS Test problem: Phillips' "famous" problem.
%
% [A,b,x] = phillips(n)
%
% Discretization of the `famous' first-kind Fredholm integral
% equation deviced by D. L. Phillips.  Define the function
%    phi(x) = | 1 + cos(x*pi/3) ,  |x| <  3 .
%             | 0               ,  |x| >= 3
% Then the kernel K, the solution f, and the right-hand side
% g are given by:
%    K(s,t) = phi(s-t) ,
%    f(t)   = phi(t) ,
%    g(s)   = (6-|s|)*(1+.5*cos(s*pi/3)) + 9/(2*pi)*sin(|s|*pi/3) .
% Both integration intervals are [-6,6].
%
% The order n must be a multiple of 4.

% Reference: D. L. Phillips, "A technique for the numerical solution
% of certain integral equations of the first kind", J. ACM 9
% (1962), 84-97.

% Discretized by Galerkin method with orthonormal box functions.

% Per Christian Hansen, IMM, 09/17/92.

% Check input.
if (rem(n,4)~=0), error('The order n must be a multiple of 4'), end

% Compute the matrix A.
h = 12/n; n4 = n/4; r1 = zeros(1,n);
c = cos((-1:n4)*4*pi/n);
r1(1:n4) = h + 9/(h*pi^2)*(2*c(2:n4+1) - c(1:n4) - c(3:n4+2));
r1(n4+1) = h/2 + 9/(h*pi^2)*(cos(4*pi/n)-1);
A = toeplitz(r1);

% Compute the right-hand side b.
if (nargout>1),
  b = zeros(n,1); c = pi/3;
  for i=n/2+1:n
    t1 = -6 + i*h; t2 = t1 - h;
    b(i) =   t1*(6-abs(t1)/2) ...
           + ((3-abs(t1)/2)*sin(c*t1) - 2/c*(cos(c*t1) - 1))/c ...
           - t2*(6-abs(t2)/2) ...
           - ((3-abs(t2)/2)*sin(c*t2) - 2/c*(cos(c*t2) - 1))/c;
    b(n-i+1) = b(i);
  end
  b = b/sqrt(h);
end

% Compute the solution x.
if (nargout==3)
  x = zeros(n,1);
  x(2*n4+1:3*n4) = (h + diff(sin((0:h:(3+10*eps))'*c))/c)/sqrt(h);
  x(n4+1:2*n4) = x(3*n4:-1:2*n4+1);
end

repo: REST


text:
function slsharedisp_attach(name, varargin)
%SLSHAREDISP_ATTACH Attachs to global display options
%
% $ Syntax $
%   - slsharedisp_attach(name, ...)
%
% $ Arguments $
%   - name:     the name of the invoker
%
% $ Description $
%   - slsharedisp_attach(name, ...) attachs the current function to
%     the global display. You can specify the following properties
%       - 'show':       whether to show process information
%                       (default = true)
%       - 'indent':     the relative indent (default = 0)
%     Please note that the actual properties set to the current displayer
%     is determined by both the properties set here and the properties
%     of the previous displayer in stack:
%       new.show = current.show && previous.show;
%       new.indent = current.indent + previous.indent;
%
% $ History $
%   - Created by Dahua Lin, on Aug 29, 2006
%

%% parse and verify input arguments

opts.show = true;
opts.indent = 0;
opts = slparseprops(opts, varargin{:});
opts.name = name;

%% Main

global GLOBAL_SHARE_DISPLAYER;

s = GLOBAL_SHARE_DISPLAYER;

if isempty(s)
    s = create_displayer([], opts);
else
    n = length(s);
    s(n+1, 1) = create_displayer(s(n), opts);
end

GLOBAL_SHARE_DISPLAYER = s;

%% Core function

function dp = create_displayer(prevdp, curdp)

dp.name = curdp.name;

if isempty(prevdp)
    dp.indent = curdp.indent;
    dp.show = curdp.show;
    dp.indentstep = 4;
else
    dp.indent = prevdp.indent + curdp.indent;
    dp.show = prevdp.show && curdp.show;
    dp.indentstep = prevdp.indentstep;
end

author: lmthang