text:
function p=unrankTComposition(theRank,t,n,firstElMostSig)
%%UNRANKTCOMPOSITION Return the composition of the given theRank, starting
%                   from zero, of n unlabeled items in t parts. That is, a
%                   method of putting n unlabeled balls into t labeled
%                   slots. Compositions are tuples of t integers >=1 that
%                   sum to n. Compositions are given in lexicographic order
%                   where the first element is the least significant unless
%                   otherwise specified.
%
%INPUTS: theRank The theRank (position in an ordered sequence) of the
%                desired composition starting from zero. There is a total
%                of binomial(n-1,t-1) compositions.
%              t The number of slots that can hold items, >=1.
%              n The number of items that are composed into slots, >=1;
%                n>=t
% firstElMostSig An optional parameter specifying whether the first
%                element is the most or least significant. The default if
%                omitted or an empty matrix is passed is false (the first
%                element is the least significant).
%
%OUTPUTS: p A tX1 vector holding the current composition, whose elements
%           sum to n. Each element is the number of "balls" in that slot.
%           If theRank is larger than the maximum number of compositions,
%           an empty matrix is returned. The values of q can range from 0
%           to n.
%
%This is an implementation of the relation described in Chapter 7.2.1.3 of
%[1] for mapping combinations to compositions. theRank is used to obtain a
%combination in lexicographic order (p(1) is the least significant
%element) and the elements of the combination are mapped to the equivalent
%composition. However, a special case is if the input is one-dimensional,
%in which case the only composition is for theRank=0 and the value is n.
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4, Fascicle 3:
%    Generating all Combinations and Partitions, Upper Saddle River, NJ:
%    Addison-Wesley, 2009.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(firstElMostSig))
    firstElMostSig=false;
end

n=n-1;
nCombo=n;
tCombo=t-1;

if(tCombo>0)
    c=unrankCombination(theRank,nCombo,tCombo);

    if(isempty(c))
        p=[];
        return;
    end

    %Transform the combination into a valid composition.
    p=zeros(t,1);
    p(1)=c(1)+1;
    for curIdx=2:tCombo
        p(curIdx)=c(curIdx)-c(curIdx-1);
    end
    p(t)=n-c(tCombo);
else%The t=1 case has to be handled separately.
    if(theRank==0)
        p=n+1;
    else
        p=[];
    end
end

if(firstElMostSig)
   p=flipud(p(:));
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef MimPluginLabelEditBox < GemLabelEditBox
    % MimPluginLabelEditBox. Part of the gui for the TD MIM Toolkit.
    %
    %     This class is used internally within the TD MIM Toolkit to help
    %     build the user interface.
    %
    %     MimPluginLabelEditBox is used to build an edit box which
    %     interacts with the MIM GUI
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = private)
        GuiApp
        Tool
        FixToInteger = true
    end

    methods
        function obj = MimPluginLabelEditBox(parent, tool, icon, gui_app)
            obj = obj@GemLabelEditBox(parent, tool.ButtonText, tool.ToolTip, class(tool));
            obj.GuiApp = gui_app;
            obj.Tool = tool;

            [value_instance_handle, value_property_name, limits_instance_handle, limits_property_name] = tool.GetHandleAndProperty(gui_app);
            value = value_instance_handle.(value_property_name);

            if ~isempty(limits_property_name)
                limits = limits_instance_handle.(limits_property_name);
            end

            obj.EditBoxPosition = tool.EditBoxPosition;
            obj.EditBoxWidth = tool.EditBoxWidth;

            if ~isempty(obj.EditBox)
                obj.EditBox.SetText(num2str(value, '%.6g'));
            end

            obj.AddPostSetListener(value_instance_handle, value_property_name, @obj.PropertyChangedCallback);

            if ~isempty(limits_property_name)
                obj.AddPostSetListener(limits_instance_handle, limits_property_name, @obj.PropertyLimitsChangedCallback);
            end
        end

        function enabled = UpdateToolEnabled(obj, gui_app)
            enabled = obj.Tool.IsEnabled(gui_app);
        end
    end

    methods (Access = protected)
        function EditBoxCallback(obj, hObject, arg2)
            EditBoxCallback@GemLabelEditBox(obj, hObject, arg2);

            [instance_handle, value_property_name, ~, ~] = obj.Tool.GetHandleAndProperty(obj.GuiApp);

            % Get the current textbox contents
            value = str2double(obj.EditBox.Text);

            % If it's a number then update the property value
            if ~isnan(value)
                if obj.FixToInteger
                    new_value = round(value);
                else
                    new_value = value;
                end
                instance_handle.(value_property_name) = new_value;
            end

            % Get the current property value. If we didn't update the
            % property value, this should be the old value. If we did
            % update it, the property set method may have modified the
            % value.
            adjusted_value = instance_handle.(value_property_name);

            % In either case we want to update the text box to reflect the
            % actual value of the property
            if adjusted_value ~= value
                obj.EditBox.SetText(num2str(adjusted_value, '%.6g'));
            end
        end

        function PropertyChangedCallback(obj, ~, ~, ~)
            [instance_handle, value_property_name, ~, ~] = obj.Tool.GetHandleAndProperty(obj.GuiApp);
            value = instance_handle.(value_property_name);
            obj.EditBox.SetText(num2str(value, '%.6g'));
        end

        function PropertyLimitsChangedCallback(obj, ~, ~, ~)
            [~, ~, limits_instance_handle, limits_property_name] = obj.Tool.GetHandleAndProperty(obj.GuiApp);
            limits = limits_instance_handle.(limits_property_name);
            range = limits(2) - limits(1);
            if abs(range) >= 100
                obj.FixToInteger = true;
            else
                obj.FixToInteger = false;
            end
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
classdef TestSparsePyrLKOpticalFlow
    %TestSparsePyrLKOpticalFlow

    methods (Static)
        function test_1
            im1 = 255*uint8([...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 1 1 1 0 0 0 0;...
                0 0 0 1 0 1 0 0 0 0;...
                0 0 0 1 1 1 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
            ]);
            im2 = circshift(im1, [0 1]);
            alg = cv.SparsePyrLKOpticalFlow();
            pts = alg.calc(im1, im2, {[3,3]});
        end

        function test_2
            prevImg = cv.imread(fullfile(mexopencv.root(),'test','RubberWhale1.png'), ...
                'Grayscale',true, 'ReduceScale',2);
            nextImg = cv.imread(fullfile(mexopencv.root(),'test','RubberWhale2.png'), ...
                'Grayscale',true, 'ReduceScale',2);
            prevPts = cv.goodFeaturesToTrack(prevImg, 'MaxCorners',200);
            alg = cv.SparsePyrLKOpticalFlow();
            [nextPts,status,err] = alg.calc(prevImg, nextImg, prevPts);
            validateattributes(nextPts, {'cell'}, ...
                {'vector', 'numel',numel(prevPts)});
            cellfun(@(pt) validateattributes(pt, {'numeric'}, ...
                {'vector', 'numel',2}), nextPts);
            validateattributes(status, {'uint8'}, ...
                {'vector', 'binary', 'numel',numel(nextPts)});
            validateattributes(err, {'single'}, ...
                {'vector', 'real', 'numel',numel(nextPts)});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function spqr_install (metis_path, tbb)
%SPQR_INSTALL compile and install SuiteSparseQR
%
% Example:
%   spqr_install                        % compiles using ../../metis-4.0, no TBB
%   spqr_install ('/my/metis')          % using non-default path to METIS
%   spqr_install ('no metis')           % do not use METIS at all
%   spqr_install ('../../metis-4.0',1)  % with METIS, and TBB multithreading
%
% SuiteSparseQR relies on CHOLMOD, AMD, and COLAMD, and can optionally use
% CCOLAMD, CAMD, and METIS as well.  By default, CCOLAMD, CAMD,
% and METIS are used.  METIS is assumed to be in the ../../metis-4.0 directory.
%
% See http://www-users.cs.umn.edu/~karypis/metis for a copy of METIS 4.0.1.
%
% You can only use spqr_install while in the SuiteSparseQR/MATLAB directory.
%
% Multithreading based on Intel Threading Building Blocks (TBB) is optionally
% used. It is not enabled by default since it conflicts with multithreading in
% the BLAS.
%
% See also spqr, spqr_solve, spqr_qmult.

%   Copyright 2008, Timothy A. Davis
%   http://www.cise.ufl.edu/research/sparse

if (nargin < 1)
    metis_path = '../../metis-4.0' ;
end

if (nargin < 2)
    tbb = 0 ;
end

% compile SuiteSparseQR and add to the path
spqr_make (metis_path, tbb) ;
spqr_path = pwd ;
addpath (spqr_path)

fprintf ('\nThe following path has been added.  You may wish to add it\n') ;
fprintf ('permanently, using the MATLAB pathtool command.\n') ;
fprintf ('%s\n', spqr_path) ;

fprintf ('\nTo try your new mexFunctions, try this command:\n') ;
fprintf ('spqr_demo\n') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function pass = sampleTest(f, sampleOP, tol, flag)
%SAMPLETEST   Test an evaluation of input OP against a SPHEREFUN.
%
%   SAMPLETEST(F, SAMPLEOP, TOL) evaluates both the function OP and its
%   SPHEREFUN representation F at several points in its domain. The
%   difference of these values is computed, and if this is sufficiently
%   small the test passes and returns TRUE.
%   If the difference is large, it returns FALSE.
%
%   SAMPLETEST(F, SAMPLEOP, TOL, FLAG) is the same as above if FLAG = 0.
%   If FLAG = 1 then the OP is assumed to be unvectorized.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [TODO]: Describe where we evaluate? (at low discrepancy points...)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% NOTE: Once complex-valued spherefuns are allowed, we can call
% sampleTest from separableApprox.

if ( nargin < 3 )
    % Assume op is vectorized:
    flag = 0;
end

domain = f.domain;

if ( ~flag )
    % Sample at lots of points if the op is vectorized.
    n = 100;
    [xeval, yeval] = halton( n, domain );

    % Evaluate the op:
    vOp = feval(sampleOP, xeval, yeval);
else
    % sample on less points if the op is unvectorized.
    n = 20;
    [xeval, yeval] = halton( n, domain );

    % Evaluate the op:
    vOp = zeros(n , 1 );
    for jj = 1:numel(xeval)
        vOp(jj) = feval(sampleOP, xeval(jj), yeval(jj));
    end
end

%for now, set to real-valued only to match what constructor does
vOp = real(vOp);

% Evaluate the SPHEREFUN:
vFun = feval(f, xeval, yeval);

% If the TECHS evaluation differs from the op evaluation, SAMPLETEST failed:
if ( any(max(abs(vOp - vFun)) > 100*tol) )
    pass = false; % :(
else
    pass = true;  % :)
end

end

function [x, y] = halton( numpts, domain )
% Halton sequences are sequences used to generate points in space, which
% are deterministic and of low discrepancy. They appear to be random for
% many purposes.
%
% Adapted from Grady Wright's code 22nd May 2014.

% generate Halton sequences on [0,1]^2:
ndims = 2;
p = [2 3 5 7 11 13];
H = zeros(numpts, ndims);
for k = 1:ndims
    N = p(k); v1 = 0; v2 = 0:N-1; lv1 = 1;
    while ( lv1 <= numpts )
        v2 = v2(1:max(2,min(N,ceil((numpts+1)/lv1))))/N;
        [x1,x2] = meshgrid(v2,v1);
        v1 = x1+x2;
        v1 = v1(:);
        lv1 = length(v1);
    end
    H(:,k) = v1(2:numpts+1);
end
% scale [0,1]^2 to the domain of the separableApprox.
x = H(:,1);
x = (domain(2) - domain(1))*x + domain(1);
y = H(:,2);
y = (domain(4) - domain(3))*y + domain(3);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%otyp==0->cont, otyp==1->disc
function [A N]=sys_bd_dcm(llh, vel_b, Cbn, gyro, otyp, dt)
% position   (1-3)
% velocity   (4-6)
% attitude   (7-9)

%system disturbance coefs
N=zeros(9,6);
N(7:9,4:6)=-Cbn; %attitude
N(4:6,1:3)=eye(3); %velocity
N(4:6,4:6)=skew(vel_b); %velocity

%system matrix
A=zeros(9);
[Rn, Re, g, sL, cL, WIE_E]=geoparam_v000(llh);
tL=sL/cL;
Rn_h=Rn+llh(3);
Re_h=Re+llh(3);
R=diag([1/Rn_h,1/Re_h,-1]);
O=[0 1/Re_h 0;-1/Rn_h 0 0;0 -tL/Re_h 0];
vel_n=Cbn*vel_b;
wen_n=[vel_n(2)/Re_h; -vel_n(1)/Rn_h; -vel_n(2)*tL/Re_h];
wie_n=[WIE_E*cL; 0; -WIE_E*sL];
win_n=wen_n+wie_n;

%Position
A(1:3,4:6)=R*Cbn;
A(1:3,7:9)=R*skew(vel_n);
A(1,3)=-vel_n(1)/Rn_h^2;
A(2,1)=vel_n(2)*tL/Re_h/cL;
A(2,3)=-vel_n(2)/Re_h^2/cL;

%Velocity
A(4:6,4:6)=-skew(Cbn'*wie_n+gyro);
A(4:6,7:9)=Cbn'*skew([0;0;-g])+skew(vel_b)*Cbn'*skew(wie_n);
A(4:6,1)=skew(vel_b)*Cbn'*[-WIE_E*sL;0;-WIE_E*cL];

%Attitude
A(7,1)=-WIE_E*sL;
A(7,3)=-vel_n(2)/Re_h^2;
A(8,3)=vel_n(1)/Rn_h^2;
A(9,1)=-(WIE_E*cL+vel_n(2)/Re_h/cL/cL);
A(9,3)=vel_n(2)*tL/Re_h^2;
A(7:9,4:6)=O*Cbn;
A(7:9,7:9)=-skew(win_n)+O*skew(vel_n);

%discretize A;
if (otyp==0)
    A=A;
elseif(otyp==1)
    A=expm(A*dt);
end

repo: ekfmonoslam


text:
function [] = bandfiltering(z_grid,x_res,y_res,spatial_bands,save_path,ifg_based_correction,n_degree_butterworth,norm_filter_flag)
% function that computes the band pass filtered data usigng 2D FFT and a
% butterwurth function. Make sure the input data is on a regular grid and
% has no NaN values. There is no output passed in this function. Instead
% the bandfiltered data is saved for each dataset individually, e.g.
% 'dataset_1.mat' for the first dataset and increasing for the others.
% input:
% z_grid                    The data specified as a grid, i.e. a matrix.
%                           Additonal datasets can be specified by
%                           increasing the third dimention.
% x_res                     The resolution in x-direction in m of the grid
% y_res                     The resolution in y-direction in m of the grid
% spatial_bands             The spatial band that need to be filtered in m,
%                           specified as a 2 column matrix [lower upper].
%                           Multiple band filters can be specified by
%                           increasing the number or rows.
% Optional inputs:
% save_path                 The path were the output data will be saved. default
%                           is the current directory
% n_degree_butterworth      The degree of the butterwurth filter, by
%                           default this is set to be 3.
% norm_filter_flag          Normalisation of the butterwurth filter is done
%                           by default. This is to scope with the issue
%                           when the extremes of the bandfitler are too
%                           close to eachother causing only a partial
%                           amplitude passing in the selected band.
%
%     Copyright (C) 2015  Bekaert David - University of Leeds
%     Email: eedpsb@leeds.ac.uk or davidbekaert.com
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
% By Bekaert David - December 2012
% modifications:
% DB	02/2013		Include the save_path variable
% DB    03/2013     Include warnings for edge effects and to large spatial
%                   bandfilters.
% DB    03/2013     Include 1D filtering to cope with limitation of narrow
%                   datasets
% DB    05/2013     Include option to crop out a region

% optional inputs for checking results:
% figure properties
fontsize = 15;
plot_figures = 0;          % plot the figures for the n_plot_dataset dataset
save_fig =0;                % when 1 save figures in the figures folder in tha aps_p folder.
n_plot_dataset=1;           % plot the figures for this band dataset.
                            % Topography is first then the interferograms
check_flag = 0;             % some additional figures being generated when turn on
mirror_flag = 1;            % when 1, mirror the dataset such filtering effects
                            % are reduced. The mirror is based on the
                            % largest spatial filter wavelength.

max_perc_mirror_2D = 100;    % filtersize with respect to the dataset spatial
                            % dimension in percentage till what 2D
                            % filtering is allowed. When larger 1D
                            % filtering is performed in the larger
                            % dimension of X or Y.

warning_flag_perc = 10;     % Output a warning when more no mirroring is done,
                            % and the percentage of half the filter length
                            % with respect to the maximum dimension is mirrored.

warning_mirror_flag_perc = 50;
                            % Output a warning when more than this percentage
                            % with respect to the maximum dimension is mirrored.
                            % Mirroring is done by half the maximum filter size.

% getting the data from the parms_aps file
crop_flag = getparm_aps('crop_flag');
 % setting the function defaults
if nargin<6
    ifg_based_correction='n';
end
if nargin<7
    n_degree_butterworth = [];
end
if nargin<8
    norm_filter_flag = [];
end
if isempty(n_degree_butterworth)==1
   n_degree_butterworth=3;
end
if isempty(norm_filter_flag)==1
   norm_filter_flag=1;
end
if nargin<5 || isempty(save_path)==1
   save_path = './';
end
fprintf(['***Bandfiltering***\n'])
fprintf(['Using butterworth filter degree: ' num2str(n_degree_butterworth) '\n'])
if norm_filter_flag==1
    fprintf(['Normalise the butterworth filter. \n'])
end

norm_filter_flag=0
n_degree_butterworth = 10
% number of datasets
n_datasets = size(z_grid,3);
n_band_filters = size(spatial_bands,1);
if strcmp(ifg_based_correction,'y')
    h_ifg_number = n_datasets/2;
else
    h_ifg_number=1;
end

if size(spatial_bands,1)==1 && size(spatial_bands,2)==2 && spatial_bands(1,1)==0 && spatial_bands(1,2)==inf
    % this is a band filter of the whole image
        for k=1:n_datasets

            if k<=h_ifg_number && h_ifg_number~=1
                % these are interferograms
                save_name = ['bandfilter_regular_hgt_ifg_' num2str(k) '.mat'];
            elseif k<=h_ifg_number && h_ifg_number==1
                % these are the heights
                save_name = 'bandfilter_regular_hgt.mat';
            else
                % thse are interferograms
                 save_name = ['bandfilter_regular_ifg_' num2str(k-h_ifg_number) '.mat'];
            end
            data_band_out(:,:) = z_grid(:,:,k);
            clear data_band

            dimension_filter=NaN;
            % saving the bandfiltered data for eacht dataset seperately
            save([save_path, filesep, save_name],'data_band_out','x_res','y_res','spatial_bands','dimension_filter')
            clear data_band_out
        end
else

    % Starting the bandfiltering code
    if plot_figures==1
        if n_plot_dataset<=h_ifg_number && h_ifg_number~=1
            % thse are interferograms
            save_folder_str=['aps_p' filesep 'fig_bandfilter_hgt_ifg' num2str(n_plot_dataset)];
        elseif n_plot_dataset<=h_ifg_number && h_ifg_number==1
            % these are the heights
            save_folder_str=['aps_p' filesep 'fig_bandfilter_hgt'];
        else
            % these are interferograms
            save_folder_str=['aps_p' filesep 'fig_bandfilter_ifg_' num2str(n_plot_dataset-h_ifg_number)];
        end

    %     if n_plot_dataset==1
    %         save_folder_str=['aps_p' filesep 'fig_bandfilter_hgt'];
    %     else
    %         save_folder_str=['aps_p' filesep 'fig_bandfilter_ifg_' num2str(n_plot_dataset-1)];
    %     end
        if exist(save_folder_str,'dir')~=7
           mkdir(save_folder_str)
        end
    end

    % original gridsize
    data_rows_or = size(z_grid,1);
    data_columns_or = size(z_grid,2);

    % investigate in the percentage of padding that is being performed
    % below the percentage of the maximum band filter with respect to the dimension is computed
    x_samples_overlap = ceil(max(spatial_bands,[],2)./x_res/2);
    x_perc_overlap = x_samples_overlap./data_columns_or*100;
    y_samples_overlap = ceil(max(spatial_bands,[],2)./y_res/2);
    y_perc_overlap = y_samples_overlap./data_rows_or*100;

    % output a warning in case the padding is to big or when the the filter
    % effect hit the warning percentage as set at the start of the code.
    if mirror_flag==1
       warning_perc = warning_mirror_flag_perc ;
    else
       warning_perc = warning_flag_perc ;
    end
    ix_x = find(x_perc_overlap>warning_perc);
    ix_y = find(y_perc_overlap>warning_perc);
    ix = intersect(ix_x,ix_y);
    clear ix_x ix_y
    if isempty(ix)~=1
        new_linestr = repmat('\n',length(ix),1);
        outputstr = [num2str(spatial_bands(ix,:)) new_linestr];
        if mirror_flag==1
            fprintf(['***Warning: The following band filters are mirrored about ',num2str(warning_perc),' perc of the maximum dimension: \n'])
            for kk=1:size(outputstr,1)
                fprintf(outputstr(kk,:)) ;
            end
            clear outputstr kk
            fprintf(['This migth introduce arctifacts. \n'])
            fprintf(['It is recommended to limit to smaller spatial bandwidths. \n'])

        else
            fprintf(['***Warning: The following band filters are above ',num2str(warning_perc),' perc of the maximum dimension: \n'])
            for kk=1:size(outputstr,1)
                fprintf(outputstr(kk,:)) ;
            end
            clear outputstr kk
            fprintf(['Likely edge effects are introduced. \n'])
            fprintf(['Turn the mirror flag on and/or limit to smaller spatial bandwidths. \n\n'])
        end
    end
    clear ix

    % output information on which bandfilters are replaced by a 1D filter to
    % reduce edge effects from a limiting dimension
    % x/y_perc_overlap represents the maximum bandfilter size with respect to
    % the x/y dimension given as a percentage
    ix_x = find(x_perc_overlap>max_perc_mirror_2D);
    ix_y = find(y_perc_overlap>max_perc_mirror_2D);

    if isempty(ix_x)~=1 || isempty(ix_y)~=1
        x_n_cases = length(ix_x);
        y_n_cases = length(ix_y);
        if x_n_cases>y_n_cases
             fprintf(['X-direction appears limited for band filtering. \n'])
             fprintf(['1D (Y-direction) band filtering is performed for the following bands: \n'])
             new_linestr = repmat('\n',length(ix_x),1);
             outputstr = [num2str(spatial_bands(ix_x,:)) new_linestr];
             for kk=1:size(outputstr,1)
                 fprintf(outputstr(kk,:)) ;
             end

             % setting the variables for the 1D filtering
             filter_1D_Y = 1;
             filter_1D_X = 0;
             ix_1D_filter = ix_x;
             clear outputstr new_linestr kk ix_y ix_x

        elseif y_n_cases>x_n_cases
             fprintf(['Y-direction appears limited for band filtering. \n'])
             fprintf(['1D (X-direction) band filtering is performed for the following bands: \n'])
             new_linestr = repmat('\n',length(ix_y),1);
             outputstr = [num2str(spatial_bands(ix_y,:)) new_linestr];
             for kk=1:size(outputstr,1)
                 fprintf(outputstr(kk,:)) ;
             end

             % setting the variables for the 1D filtering
             filter_1D_Y = 0;
             filter_1D_X = 1;
             ix_1D_filter = ix_y;
             clear outputstr new_linestr kk ix_x ix_y

        else
             fprintf(['Both dataset appears as limited for band filtering. \n'])
             fprintf(['2D band filtering is performed, but check the following bands: \n'])
             new_linestr = repmat('\n',length(ix_y),1);
             outputstr = [num2str(spatial_bands(ix_y,:)) new_linestr];
             for kk=1:size(outputstr,1)
                 fprintf(outputstr(kk,:)) ;
             end

             % setting the variables for the 1D filtering
             filter_1D_Y = 0;
             filter_1D_X = 0;
             ix_1D_filter = [];
             clear outputstr new_linestr kk ix_x ix_y
        end
    else
         filter_1D_Y = 0;
         filter_1D_X = 0;
         ix_1D_filter = [];
    end

    if plot_figures==1
        data_original_image = z_grid(:,:,n_plot_dataset);
    end

    % mirror the edges of the grid to reduce filter effects at the edges
    % get the maximum size of the spatial filters
    if mirror_flag == 1
        fprintf('Perform mirroring to reduce filtering effects on edges \n')
        fprintf('By half the maximum filter length at each edge. \n')

        % maximum spatial wavelength
        max_band = max(max(spatial_bands));
        if max_band == inf
           max_band = 100000;        % extend the grid to a maximum of 100 km in case one goes to infinite
        end
        n_mirror_x = ceil(max_band./x_res*1.5);
        n_mirror_y = ceil(max_band./y_res*1.5);

        for k=1:n_datasets
            % padding the grid symmetric
            z_grid_new(:,:,k) = padarray(z_grid(:,:,k),[n_mirror_y n_mirror_x],'symmetric');
            % plotting an intermediate figure when requested
            if plot_figures==1 && k==n_plot_dataset
                h1= figure('name','Data symmetric padded for largest filter');
                imagesc(z_grid_new(:,:,k))
                axis equal
                axis tight

                % saving of the figure when requested
                if save_fig==1
                    fig_save_name = [save_folder_str filesep 'original_data_mirrored.eps'];
                    set(h1,'PaperPositionMode','auto')
                    print(h1,'-depsc','-r150',fig_save_name)
                    clear fig_save_name
                    close(h1)
                end
                clear h1
            end
        end
        clear z_grid
        z_grid = z_grid_new;
    else
       fprintf('No mirroring performed. Filter artifacts will become more persistent for larger wavelengths! \n')
    end

    % size of the dataset
    data_rows = size(z_grid,1);
    data_columns = size(z_grid,2);

    % Sampling frequency follows from the resolution
    fs_rows = 1/y_res;          % rows sampling frequency [1/m]
    fs_columns = 1/x_res;   	% columns sampling frequency [1/m]

    % rows and columns of the data such they are a power of 2
    % this will speed up the fft and will autmoatically padd
    % the data matrix with zeros where needed
    data_rows_new = 2.^nextpow2(data_rows);
    data_columns_new = 2.^nextpow2(data_columns);

    if plot_figures==1
        if mirror_flag==1
            % computation of the axis extremes based on the resolution given
            X_lims_fig = [0 data_columns_or*x_res];     % axis limits in [m]
            Y_lims_fig = [0 data_rows_or*y_res];        % axis limits in [m]
        else
            X_lims_fig = [0 data_columns*x_res];        % axis limits in [m]
            Y_lims_fig = [0 data_rows*y_res];           % axis limits in [m]
        end
        % computation of the axis extremes based on the resolution given
        % Asummed to be symmetric padded. In case no padding is done, this will
        % be equal to the orginal dataset
        X_lims_fig_syix_1D_filterm = [0 data_columns*x_res];        % axis limits in [m]
        Y_lims_fig_sym = [0 data_rows*y_res];           % axis limits in [m]
    end

    % Loop over each dataset. Invert to the freq domain once and reinvert for
    % each bandfilter. i.e. There will be a loop over the different
    % bandfilters. For each dataset the data is being saved in a matfile.
    for k=1:n_datasets
        if k==n_plot_dataset && plot_figures==1
            plot_fig_flag = 1;
        else
            plot_fig_flag = 0;
        end

        fprintf(['Progress: ' num2str(k) '/' num2str(n_datasets) ,' done \n'])
        % when chosen visualise the dataset
        if plot_fig_flag==1
            h1 = figure('name','Original image');
            imagesc(X_lims_fig,Y_lims_fig,data_original_image)
            cc = colorbar;
            axis xy
            axis equal
            axis tight
            set(gca,'fontsize',fontsize)
            xlabel('Distance [m]','fontsize',fontsize)
            ylabel('Distance [m]','fontsize',fontsize)
            title('Original image','fontsize',fontsize)
            colorlimits = get(cc,'YLim');

            % saving of the figure when requested
            if save_fig==1
                fig_save_name = [save_folder_str filesep 'input.eps'];
                set(h1,'PaperPositionMode','auto')
                print(h1,'-depsc','-r150',fig_save_name)
                clear fig_save_name
                close(h1)
            end
            clear h1

        end

        % converting to the frequency domain by FFT
        data_freq_complex = fft2(z_grid(:,:,k),data_rows_new,data_columns_new);

        % shifting the spectrum around zero freq
        % matlab shifts with 1 pixel off for even number columns and rows
        data_freq_complex = fftshift(data_freq_complex);
        data_freq_complex_new = [data_freq_complex  data_freq_complex(:,1)];
        data_freq_complex_new = [data_freq_complex_new ; data_freq_complex_new(1,:)];
        clear data_freq_complex

        % Frequency domain: sampling resolution follows from the Nyquist requency (spampling freq/samples)
        % freq resolution of the bins in the freq domain
        fs_rows = 1/data_rows_new*1/y_res;
        fs_columns = 1/data_columns_new*1/x_res;

        % Computing the frequencies for the axis in the freq domain figures.
        freq_rows_fig_vector = [-(data_rows_new-1)/2-0.5:1:(data_rows_new-1)/2+0.5].*fs_rows;
        freq_columns_fig_vector= [-(data_columns_new-1)/2-0.5:1:(data_columns_new-1)/2+0.5].*fs_columns;

        if plot_fig_flag==1  && check_flag==1
            figure('name','Spectrum centralised around zero freq')
            imagesc(freq_columns_fig_vector,freq_rows_fig_vector,abs(data_freq_complex_new))
            xlabel('Xfreq [1/m]','fontsize',fontsize)
            ylabel('Yfreq [1/m]','fontsize',fontsize)
            title('Spectrum centralised around zero freq','fontsize',fontsize)
            colorbar
            set(gca,'fontsize',fontsize)
            axis equal
            axis tight

            figure('name','Phase centralised around zero freq')
            imagesc(freq_columns_fig_vector,freq_rows_fig_vector,angle(data_freq_complex_new))
            xlabel('Xfreq [1/m]','fontsize',fontsize)
            ylabel('Yfreq [1/m]','fontsize',fontsize)
            title('Phase centralised around zero freq','fontsize',fontsize)
            colorbar
            set(gca,'fontsize',fontsize)
            axis equal
            axis tight
        end

        %% Computing the grid of frequencies
        freq_rows_fig_matrix = repmat(freq_rows_fig_vector',1,data_columns_new+1);
        freq_columns_fig_matrix = repmat(freq_columns_fig_vector,data_rows_new+1,1);
        % Set the frequencies for 1D or 2D filtering
        if filter_1D_X==1
            % 1D X-direction filter
            freq_fig_matrix_1D = freq_columns_fig_matrix;
            if plot_fig_flag==1 && check_flag==1
                figure('name','Frequency grid 1D X-direction filter')
                imagesc(freq_fig_matrix_1D)
                xlabel('X','fontsize',fontsize)
                ylabel('Y','fontsize',fontsize)
                cc = colorbar;
                set(gca,'fontsize',fontsize)
                xlabel(cc,'Freq [1/m]','fontsize',fontsize)
            end
        elseif filter_1D_Y==1
            % 1D Y-direction filter
            freq_fig_matrix_1D = freq_rows_fig_matrix;
            if plot_fig_flag==1 && check_flag==1
                figure('name','Frequency grid 1D Y-direction filter')
                imagesc(freq_fig_matrix_1D)
                xlabel('X','fontsize',fontsize)
                ylabel('Y','fontsize',fontsize)
                cc = colorbar;
                set(gca,'fontsize',fontsize)
                xlabel(cc,'Freq [1/m]','fontsize',fontsize)
            end
        end
        % 2D filter
        freq_fig_matrix_2D = sqrt(freq_rows_fig_matrix.^2+freq_columns_fig_matrix.^2);
        if plot_fig_flag==1 && check_flag==1
            figure('name','Frequency grid')
            imagesc(freq_fig_matrix_2D)
            xlabel('X','fontsize',fontsize)
            ylabel('Y','fontsize',fontsize)
            cc = colorbar;
            set(gca,'fontsize',fontsize)
            xlabel(cc,'Freq [1/m]','fontsize',fontsize)
        end

        %% Loop over the different bandfilters
        % initialisation of the output data_band variable
        % This is after the symmetric padding has been removed again.
        data_band_out = NaN([data_rows_or data_columns_or n_band_filters]);
        dimension_filter = NaN([n_band_filters 1]);
        for kk=1:n_band_filters
            % output to the screen
            fprintf(['Spatial bandfilter: ',num2str(spatial_bands(kk,1)), '\t - ',num2str(spatial_bands(kk,2)),' \t m \t'])

            % convert spatial wavelength band to spatial frequency band
            f_band = 1./sort(spatial_bands(kk,:));	% frequency band [1/m]

            % Checking if the band filter is set as 1D or 2D filter
            filter_1D = find(kk==ix_1D_filter);
            if isempty(filter_1D)~=1
                % do 1D filtering
                % computing the function of the bandfilter
                H_low = 1./(1+(freq_fig_matrix_1D./f_band(1)).^(2*n_degree_butterworth));
                H_high = 1./(1+(freq_fig_matrix_1D./f_band(2)).^(2*n_degree_butterworth));
                dimension_filter(kk) = 1;
                fprintf(['(1D filter)\n'])
            else
                % do 2D filtering
                % computing the function of the bandfilter
                H_low = 1./(1+(freq_fig_matrix_2D./f_band(1)).^(2*n_degree_butterworth));
                H_high = 1./(1+(freq_fig_matrix_2D./f_band(2)).^(2*n_degree_butterworth));
                dimension_filter(kk) = 2;
                fprintf(['(2D filter)\n'])
            end

            % Correct for the case of a filter having inf in its range
            ix_center_spectrum = (size(H_high)-1)./2+1;
            if isnan(H_high(ix_center_spectrum(1),ix_center_spectrum(2)))
                H_high(ix_center_spectrum(1),ix_center_spectrum(2))=0;
            end
            if isnan(H_low(ix_center_spectrum(1),ix_center_spectrum(2)))
                H_low(ix_center_spectrum(1),ix_center_spectrum(2))=0;
            end
            % Computing the band pass filter
            H_Butterworth =  (H_low - H_high);
            H_Butterworth_norm = H_Butterworth./max(max(H_Butterworth));

            % plotting when requested
            if plot_fig_flag==1
                % plotting the filter around zero for the rows
                h1 = figure('name','Rows freq filter');
                subplot(4,1,1)
                plot(freq_rows_fig_vector,H_low(:,ix_center_spectrum(2)),'b.-')
                ylim([0 1])
                title(['freq : ' num2str(f_band(1)) ])
                subplot(4,1,2)
                plot(freq_rows_fig_vector,H_high(:,ix_center_spectrum(2)),'b.-')
                ylim([0 1])
                title(['freq : ' num2str(f_band(2)) ])
                subplot(4,1,3)
                plot(freq_rows_fig_vector,H_Butterworth(:,ix_center_spectrum(2)),'b.-')
                ylim([0 1])
                title(['Band freq' ])
                subplot(4,1,4)
                plot(freq_rows_fig_vector,H_Butterworth_norm(:,ix_center_spectrum(2)),'b.-')
                ylim([0 1])
                title(['Norm band freq' ])

                % saving of the figure when requested
                if save_fig==1
                    fig_save_name = [save_folder_str filesep 'bandfilter_' num2str(spatial_bands(kk,1)) '_'  num2str(spatial_bands(kk,2)) 'm.eps'];
                    set(h1,'PaperPositionMode','auto')
                    print(h1,'-depsc','-r150',fig_save_name)
                    clear fig_save_name
                    close(h1)
                end
                clear h1

            end

            if norm_filter_flag ==1
                H_Butterworth = H_Butterworth_norm;
            end

            % multiplying filter in frequency domain.
            % data_freq_complex_new = real(data_freq_complex_new).*H_Butterworth+1i.*angle(data_freq_complex_new);
            data_freq_complex_new_band = data_freq_complex_new.*H_Butterworth;
            if plot_fig_flag==1  && check_flag==1
                figure('name','New spectrum centralised around zero freq')
                imagesc(freq_columns_fig_vector,freq_rows_fig_vector,abs(data_freq_complex_new_band))
                xlabel('Xfreq [1/m]','fontsize',fontsize)
                ylabel('Yfreq [1/m]','fontsize',fontsize)
                title('Spectrum centralised around zero freq','fontsize',fontsize)
                colorbar
                set(gca,'fontsize',fontsize)
                axis equal
                axis tight

                figure('name','New phase centralised around zero freq')
                imagesc(freq_columns_fig_vector,freq_rows_fig_vector,angle(data_freq_complex_new_band))
                xlabel('Xfreq [1/m]','fontsize',fontsize)
                ylabel('Yfreq [1/m]','fontsize',fontsize)
                title('Phase centralised around zero freq','fontsize',fontsize)
                colorbar
                set(gca,'fontsize',fontsize)
                axis equal
                axis tight
            end

            % Shifting the spectrum such the center frequency is at the top left corner
            data_freq_complex_new_band(:,end) = [];
            data_freq_complex_new_band(end,:) = [];
            data_freq_complex_new_band = ifftshift(data_freq_complex_new_band);

            % Inverting back to the spatial domain
            data_band = ifft2(data_freq_complex_new_band,data_rows_new,data_columns_new);

            % Removing the padded region introduced when putting the number of lines and rows 2^n
            if data_rows-data_rows_new<0			% for the rows
                data_band(data_rows+1:end,:)=[];
            end
            if data_columns-data_columns_new<0		% for the columns
                data_band(:,data_columns+1:end)=[];
            end

            % Keeping only the real part
            data_band = real(data_band);

            % removing any mirrored image from the data
            if mirror_flag == 1
                data_band_or = data_band;
                data_band([1:n_mirror_y],:)=[];
                data_band([end-n_mirror_y+1:end],:)=[];
                data_band(:,[1:n_mirror_x])=[];
                data_band(:,[end-n_mirror_x+1:end])=[];
            end

            % plotting the final results when requested
            if plot_fig_flag==1
                h1= figure('name','output image');
                imagesc(X_lims_fig,Y_lims_fig,data_band)
                % limit the colorbar to 95% bounds
                data_sorted = sort(reshape(data_band,[],1));
                ix = round([0.025 0.975].*length(data_sorted));
                colorlimits = [data_sorted(ix)];
                caxis(colorlimits)
                colorbar
                set(gca,'fontsize',fontsize)
                xlabel('Distance [m]','fontsize',fontsize)
                ylabel('Distance [m]','fontsize',fontsize)
                title(['Band filter: ' , num2str(spatial_bands(kk,1)) ,' -- ' ,num2str(spatial_bands(kk,2)),' m'],'fontsize',fontsize)
                axis equal
                axis tight
                axis xy

                % saving of the figure when requested
                if save_fig==1
                    fig_save_name = [save_folder_str filesep 'output_' num2str(spatial_bands(kk,1)) '_'  num2str(spatial_bands(kk,2)) 'm.eps'];
                    set(h1,'PaperPositionMode','auto')
                    print(h1,'-depsc','-r150',fig_save_name)
                    clear fig_save_name
                    close(h1)
                end
                clear h1

            end

            % Storing the output data
            data_band_out(:,:,kk) = data_band;
            clear data_band
        end

        if k<=h_ifg_number && h_ifg_number~=1
            % thse are interferograms
            save_name = ['bandfilter_regular_hgt_ifg_' num2str(k) '.mat'];
        elseif k<=h_ifg_number && h_ifg_number==1
            % these are the heights
            save_name = 'bandfilter_regular_hgt.mat';
        else
            % thse are interferograms
             save_name = ['bandfilter_regular_ifg_' num2str(k-h_ifg_number) '.mat'];
        end

        % saving the bandfiltered data for eacht dataset seperately
        save([save_path, filesep, save_name],'data_band_out','x_res','y_res','spatial_bands','n_degree_butterworth','norm_filter_flag','dimension_filter')
        clear data_band_out
    end
end

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
classdef dat_wc_reader < handle
	properties
        sr
        min_index
        max_segments
        raw_filename
        sr_infile
        t0_segments
        max_index
        segmentLength
        spikes_file
    end
	methods
        function obj = dat_wc_reader(par, raw_filename)
            obj.sr = [];
            obj.max_segments = [];
            obj.raw_filename = raw_filename;
            obj.spikes_file = false;
            data = csvread(raw_filename);

            obj.sr_infile = false;
            if isfield(par,'sr')
                obj.sr = par.sr;
            end
            obj.min_index = floor(par.tmin * obj.sr);                   %min time to read (in micro-sec)

            if strcmp(par.tmax,'all')
                obj.max_index = max(size(data));
            else
                obj.max_index = ceil(par.tmax * obj.sr);
            end
            t0 = (obj.min_index-1)/obj.sr*1000;
            n = obj.max_index -  obj.min_index;

            %Segments the data in par.segments pieces
            obj.max_segments = ceil(n/ obj.sr / ...
                    (par.segments_length * 60));         %number of segments in which data is cutted

            obj.segmentLength = floor (n/obj.max_segments);

            obj.t0_segments = ones(1,obj.max_segments);
            obj.t0_segments(1) = t0;
            for i = 2:obj.max_segments
                obj.t0_segments(i) = obj.t0_segments(i-1) + obj.segmentLength/obj.sr*1000;
            end

        end

        function [sr,max_segments,with_raw,with_spikes] = get_info(obj)
            if isempty(obj.max_segments)
                with_raw = false;
                max_segments = 0;
            else
                with_raw = true;
                max_segments = obj.max_segments;
            end

            with_spikes = obj.spikes_file;
            if obj.sr_infile
                sr = obj.sr;
            else
                sr = [];
            end
        end

        function index_ts = index2ts(obj,index,i)
            index_ts = (index)/obj.sr*1000 + obj.t0_segments(i);
        end

        function x = get_segment(obj,i)

            if i ~= obj.max_segments
                lastsample = obj.min_index+obj.segmentLength*i;
            else
                lastsample = obj.max_index;
            end
            firstsample = obj.min_index+obj.segmentLength*(i-1)+1;

            x = csvread(obj.raw_filename); %it's sub optimal but the best way to handle inconsistent .dat files
            x = x(firstsample:lastsample);
        end
    end
end

author: csn-le


text:
function mask = inplaneOverlayMask(view, img, hImg);
%
% mask = inplaneOverlayMask(view, img, <hImg>);
%
% Compute a binary mask indicating which pixels in the current
% overlay image for a view exceed the thresholds set by the UI.
%
% img is a data matrix containing the overlay data.
%
% hImg is an optional handle to an overlay image. If this is passed,
% will update the image's Alpha map such that the pixels exceeding
% threshold are opaque and everything else is transparent.
%
% ras, 05/2006.
if nargin<1, view = getCurView; end

if ~exist('img', 'var') | isempty(img),
    if checkfields(view, 'ui', 'image')
        img = view.ui.image;
    else
        error('Need an overlay data argument (img).')
    end
end

% Select pixels that satisfy cothresh, phWindow, and mapWindow
cothresh = getCothresh(view);
phWindow = getPhWindow(view);
mapWindow = getMapWindow(view);

pts = [];
% if ~isempty(img)
%     pts = ones(size(img));
%     curCo = cropCurSlice(view,'co',slice);
%     curPh = cropCurSlice(view,'ph',slice);
%     curMap = cropCurSlice(view,'map',slice);
%
%     if ~isempty(curCo) & cothresh>0
%         ptsCo = curCo > cothresh;
%         pts = pts & ptsCo;
%     end
%
%     if ~isempty(curPh)
%         if diff(phWindow) > 0
%             ptsPh = (curPh>=phWindow(1) & curPh<=phWindow(2));
%         else
%             ptsPh = (curPh>=phWindow(1) | curPh<=phWindow(2));
%         end
%         pts = pts & ptsPh;
%     end
%
%     if ~isempty(curMap)
%         ptsMap = (curMap>=mapWindow(1) & curMap<=mapWindow(2));
%         pts = pts & ptsMap;
%     end
% end

switch view.ui.displayMode
    case 'anat', mask = logical(zeros(size(img)));
    case 'map',
        mask = (img >= mapWindow(1)) & (img <= mapWindow(2));
    case 'amp',
        % tough ...
        mask = (img >= mapWindow(1)) & (img <= mapWindow(2));
    case 'ph',
        mask = (img >= mapWindow(1)) & (img <= mapWindow(2));
    case 'co',
        mask = (img >= mapWindow(1)) & (img <= mapWindow(2));
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function F = in_fread_bci2000(sFile, SamplesBounds)
% IN_FREAD_BCI2000: Read a block of recordings from a BCI2000 .dat file
%
% Uses library: https://www.bci2000.org/mediawiki/index.php/User_Reference:Matlab_MEX_Files

% @=============================================================================
% 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.
% =============================================================================@
%
% Author: Francois Tadel 2022

% Parse inputs
if (nargin < 2) || isempty(SamplesBounds)
    if ~isempty(sFile.epochs)
        SamplesBounds = round(sFile.epochs(iEpoch).times .* sFile.prop.sfreq);
    else
        SamplesBounds = round(sFile.prop.times .* sFile.prop.sfreq);
    end
end

% Install plugin BCI2000
if ~exist('load_bcidat', 'file')
    [isInstalled, errMsg] = bst_plugin('Install', 'bci2000');
    if ~isInstalled
        error(errMsg);
    end
end

% Read signals
F = load_bcidat(sFile.filename, SamplesBounds + 1, '-calibrated')';

author: brainstorm-tools


text:
function c = tapas_condhalluc_obs2_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the response model used to analyze data from conditioned
% hallucination paradigm by Powers & Corlett
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The rationale for this model is as follows:
%
% TO BE DESCRIBED...
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2016 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'tapas_condhalluc_obs2';

% Sufficient statistics of Gaussian parameter priors

% Beta
c.logbemu = log(48);
c.logbesa = 1;

% Nu
c.lognumu = log(1);
c.lognusa = 1;

% Gather prior settings in vectors
c.priormus = [
    c.logbemu,...
    c.lognumu,...
         ];

c.priorsas = [
    c.logbesa,...
    c.lognusa,...
         ];

% Model filehandle
c.obs_fun = @tapas_condhalluc_obs2;

% Handle to function that transforms observation parameters to their native space
% from the space they are estimated in
c.transp_obs_fun = @tapas_condhalluc_obs2_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
function spm_bms_display_ROI (BMS,mask,method)
% display results from BMS in a region of interest (ROI)
% FORMAT spm_bms_display_ROI (BMS,mask,method)
%
% Input:
% BMS    - BMS.mat file
% mask   - region of interest image
% method - inference method (FFX or RFX)
%__________________________________________________________________________
% Copyright (C) 2009-2011 Wellcome Trust Centre for Neuroimaging

% Maria Joao Rosa
% $Id: spm_bms_display_ROI.m 5219 2013-01-29 17:07:07Z spm $

% Find graphics window
% -------------------------------------------------------------------------
Fgraph = spm_figure('GetWin','Graphics');

% Input
% -------------------------------------------------------------------------
if nargin<1
    BMS    = spm_select(1,'^BMS.mat$','select BMS.mat files');
end
if nargin<2
    mask   = spm_select(1,'image','select ROI file');
end
if nargin<3
    method = spm_input('Inference method',+1,'b','FFX|RFX',['ffx';'rfx']);
end

mask_image = spm_vol(mask);         % Mask image Vol

% Nb. of subjects and models
% -------------------------------------------------------------------------
if isfield(BMS.map,method)
   switch method
       case 'ffx'
            data   = BMS.map.ffx.data;
       case 'rfx'
            data   = BMS.map.rfx.data;
   end
else
  msgbox(sprintf('Error: no %s analysis in current BMS.mat!',method))
  return
end

nsubjs    = size(data,2);
nmodels   = size(data{1}(1).mod_map,1);
nsess     = size(data{1},2);

% Sort out log-evidence images dimensions
% -------------------------------------------------------------------------
Vol_models(1,1) = spm_vol(data{1}.mod_map(1));

first_vol       = Vol_models(1,1);
M               = first_vol{1}.mat;
DIM             = first_vol{1}.dim(1:3)';

xdim            = DIM(1);
ydim            = DIM(2);
zdim            = DIM(3);
[xords,yords]   = ndgrid(1:xdim,1:ydim);
xords           = xords(:)';
yords           = yords(:)';
I               = 1:xdim*ydim;
zords_init      = ones(1,xdim*ydim);

% Loop through data
% -------------------------------------------------------------------------
for i = 1:nmodels,
    for s = 1:nsubjs,
        for se = 1:nsess,
             Vol_models(s,i,se) = spm_vol(data{s}(se).mod_map(i));
        end
    end
end

log_ev_roi = zeros(nsubjs,nmodels);
nvox_total = 0;

% Loop through image slices
% -------------------------------------------------------------------------
for z = 1:zdim,

    j = NaN(xdim,ydim);                         % Init. image values

    zords   = z*zords_init;                     % Slice z
    xyz     = [xords(I); yords(I); zords(I)];   % Slice coordinates
    nVox    = size(xyz,2);                      % Nb. of voxels per slice

    % Voxels inside mask
    mask_xyz  = mask_image.mat\M*[xyz(:,1:nVox);ones(1,nVox)];
    gamma     = spm_get_data(mask_image,mask_xyz);
    b         = find(gamma>0.5);            % Voxels in the mask

    z_models        = NaN(nsubjs,nmodels,nVox);       % Data
    z_models(1,1,:) = spm_get_data(first_vol{1},xyz); % Data: all subs/mods
    non_nan         = find(~isnan(z_models(1,1,:)));  % Voxels ~NaN

    % Find voxels ~NaN and sum sessions
    % ---------------------------------------------------------------------
    for s = 1:nsubjs,
        for k = 1:nmodels,
                sum_tmp_data    = [];
            for ns = 1:nsess,
                tmp_data        = Vol_models(s,k,ns);
                sum_tmp_data    = [sum_tmp_data; spm_get_data(tmp_data{1},xyz)];
            end
                z_models(s,k,:) = sum(sum_tmp_data,1);
                non_nani        = find(~isnan(z_models(s,k,:)));
                non_nan         = intersect(non_nan,non_nani);
        end
    end

    % Voxels to be analysed
    non_nan = intersect(non_nan,b);
    Nvoxels = length(non_nan);

    if Nvoxels > 0

        nvox_total = nvox_total + Nvoxels;
        % Do BMS in all voxels of slice z
        for n = 1:Nvoxels,
            log_ev_tmp = z_models(:,:,non_nan(n));
            % Group BMS
            log_ev_roi = log_ev_roi + log_ev_tmp;

        end
    end
end % Loop over slices

% Method
% -------------------------------------------------------------------------
switch method

    % Fixed Effects
    % ---------------------------------------------------------------------
    case 'ffx',

        F    = sum(log_ev_roi) - min(log_ev_roi);
        i    = F < (max(F) - 32);
        P    = F;
        P(i) = max(F) - 32;
        P    = P - min(P);
        P    = exp(P);
        P    = P/sum(P);

        if  isfield(BMS.map,'ffx')

            % Bar plot
            figure(Fgraph);
            spm_results_ui('Clear',Fgraph);

            hvox   = axes('Position',[0.25 0.15 0.5 0.25],'Parent',...
                Fgraph,'Visible','off');

            bar(1:nmodels,P)
            set(gca,'XTick',1:nmodels)
            set(gca,'XTickLabel',1:nmodels)
            set(gca,'YLim',[0 1])
            ylabel('Posterior Model Probability','Fontsize',12)
            xlabel('Models','Fontsize',12)
            title({'Fixed-effects BMS';''},...
                'Fontsize',12);
            axis square
            grid on

            return

        else

            msgbox('Error: no FFX analysis in current BMS.mat!')
            return

        end

    % Random Effects
    % ---------------------------------------------------------------------
    case 'rfx',

        nsamps           = 1e3;
        [alpha,exp_r,xp] = spm_BMS(log_ev_roi,nsamps,0,0,1);

        if  isfield(BMS.map,'rfx')

            % Bar plots
            figure(Fgraph);
            spm_results_ui('Clear',Fgraph);

            hvox   = axes('Position',[0.55 0.18 0.30 0.20],'Parent',...
                    Fgraph,'Visible','off');

            bar(1:nmodels,xp)
            set(gca,'XTick',1:nmodels)
            set(gca,'XTickLabel',1:nmodels)
            set(gca,'YLim',[0 1])
            ylabel('Exceedance Probability','Fontsize',12)
            xlabel('Models','Fontsize',12)
            title({'Random-effects BMS';''},'Fontsize',12)
            axis square
            grid on

            hvox   = axes('Position',[0.16 0.18 0.30 0.20],'Parent',...
                    Fgraph,'Visible','off');

            bar(1:nmodels,exp_r)
            set(gca,'XTick',1:nmodels)
            set(gca,'XTickLabel',1:nmodels)
            set(gca,'YLim',[0 1])
            ylabel('Expected Posterior Probability','Fontsize',12)
            xlabel('Models','Fontsize',12)
            title({'Random-effects BMS';''},'Fontsize',12)
            axis square
            grid on

            return

        else

            msgbox('Error: no RFX analysis in current BMS.mat!')
            return

        end

end

author: spm


text:
% % This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% % Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% % Copyright: ETH Zurich, 2014
% % Author: Romano Cicchetti
%
% function create_configuration()
%
%     algorithm = 'weiss';
%     configuration = 'initial';
%
%     % create a new default configuration
%
%     % load default values in yaml file
%     yaml_file = strcat('input/default_values/', algorithm, '.yaml');
%     default = ReadYaml(yaml_file);
%
%     % set default values
%     field_names = fieldnames(default);
%     for field_name = field_names'
%         default.(field_name{1}) = cell2mat(default.(field_name{1}));
%     end
%
%     % save default values in matlab file and a new yaml file
%     path_to_configuration = strcat('input/autogen/configurations/', algorithm, '_', configuration);
%     if ~exist(path_to_configuration, 'dir')
%         mkdir(path_to_configuration);
%     end
%     save(strcat(path_to_configuration, '/default_values.mat'), 'default');
%     copyfile(yaml_file, strcat(path_to_configuration, '/default_values.yaml'));
%
% end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function im_merged = mergeAvg(varargin)
% Average merging on channels
% Support a changing number of inputs
% Merge via Euclidean norm
im_merged = sqrt(sum(cat(3, varargin{:}).^2, 3));
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
function obj_moveBack(ho)
%OBJ_MOVEBACK - Move graphic object to the background
%
%Synposis:
% obj_moveBack(<HO>)
%
%Input:
% HO: handle of graphic object, default: last plotted object

if ~exist('ho', 'var'),
  ho= gco;
elseif length(ho)>1,
  for h= ho(:)',
    obj_moveBack(h);
  end
  return;
end

hp= get(ho, 'parent');
hc= get(hp, 'children');
if ~exist('ho', 'var') || isempty(ho),
  hi= length(hc);
else
  hi= find(hc==ho);
  if isempty(hi),
    error('handle not found in current axes');
  end
end

set(hp, 'children', [hc([1:hi-1, hi+1:end, hi])]);

author: PatternRecognition

repo: OpenBMI


text:
function sudoku_test02 ( )

%*****************************************************************************80
%
%% SUDOKU_TEST02 demonstrates SUDOKU_DIGIT_PERMUTE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SUDOKU_TEST02:\n' );
  fprintf ( 1, '  SUDOKU_DIGIT_PERMUTE permutes the digits\n' );
  fprintf ( 1, '  of a Sudoku puzzle.\n' );

  s = [ ...
       0,0,1,9,0,0,0,0,8; ...
       6,0,0,0,8,5,0,3,0; ...
       0,0,7,0,6,0,1,0,0; ...
       0,3,4,0,9,0,0,0,0; ...
       0,0,0,5,0,4,0,0,0; ...
       0,0,0,0,1,0,4,2,0; ...
       0,0,5,0,7,0,9,0,0; ...
       0,1,0,8,4,0,0,0,7; ...
       7,0,0,0,0,9,2,0,0];

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here is the initial Sudoku puzzle:\n' );
  fprintf ( 1, '\n' );

  sudoku_print ( s );

  permutation = [ 2, 3, 4, 5, 6, 7, 8, 9, 1 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here is the permutation:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Original  Permuted\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 9
    fprintf ( 1, '  %8d  %8d\n', i, permutation(i) );
  end

  s = sudoku_digit_permute ( permutation, s );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here is the permuted Sudoku:\n' );
  fprintf ( 1, '\n' );

  sudoku_print ( s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef Server < handle
  properties(SetAccess=private, GetAccess=private)
    heightmaps=containers.Map('KeyType', 'uint32', 'ValueType', 'any');
    map_id_last_used=[];
    map_id_last_added=[];
    max_stored_maps=5;
  end

  methods
    function obj = Server()
    end

    function region = getCSpaceRegionAtIndex(obj, i0, yaw, collision_model, varargin)
      p = inputParser();
      p.addRequired('i0', @isnumeric);
      p.addRequired('yaw', @isnumeric);
      p.addRequired('collision_model', @(x) isa(x, 'iris.terrain_grid.CollisionModel'));
      p.addParamValue('map_id', -1, @isnumeric);
      p.addParamValue('xy_bounds', iris.Polyhedron(zeros(0,2),zeros(0,1)), @(x) isa(x, 'iris.Polyhedron'));
      p.addParamValue('plane_distance_tolerance', 0.025, @isnumeric);
      p.addParamValue('plane_angle_tolerance', 10 * pi/180, @isnumeric);
      p.addParamValue('excluded_grid', []);
      p.addParamValue('debug', false);
      p.addParamValue('error_on_infeasible_start', true);
      p.parse(i0, yaw, collision_model, varargin{:});
      options = p.Results;

      if options.map_id == -1
        heightmap = obj.getHeightmap(obj.map_id_last_added(end));
      else
        heightmap = obj.getHeightmap(options.map_id);
      end

      x0 = heightmap.X(i0);
      y0 = heightmap.Y(i0);
      z0 = heightmap.Z(i0);
      p0 = [x0;y0;z0];
      n0 = heightmap.normals(:,i0);

      sz = size(heightmap.X);

      dist_to_plane = abs(n0'*[reshape(heightmap.X, 1, []); reshape(heightmap.Y, 1, []); reshape(heightmap.Z, 1, [])] - n0'*p0);
      dist_to_plane = reshape(dist_to_plane, sz);
      dist_mask = abs(dist_to_plane < options.plane_distance_tolerance);

      normal_product = n0' * heightmap.normals;
      normal_product = reshape(normal_product, sz);
      normal_angle_mask = normal_product > cos(options.plane_angle_tolerance);

      plane_mask = dist_mask & normal_angle_mask;
      if ~isempty(options.excluded_grid)
        if ~options.excluded_grid(i0)
          error('IRIS:TerrainSegmentation:BadSeed', 'the seed point is marked as infeasible in the excluded_grid option');
        end
        plane_mask = plane_mask & options.excluded_grid;
      end

      if ~plane_mask(i0)
        error('IRIS:TerrainSegmentation:SelectedPointInfeasible', 'cannot create a region around this point.');
      end

      boundary_mask = logical(iris.terrain_grid.component_boundary(plane_mask, i0));

      if options.debug
        figure(4)
        clf
        subplot(321)
        imshow(dist_mask, 'InitialMagnification', 'fit');
        title('dist')
        subplot(323)
        imshow(normal_angle_mask, 'InitialMagnification', 'fit');
        title('normal');
        subplot(325)
        imshow(options.excluded_grid, 'InitialMagnification', 'fit');
        title('existing');
        subplot(322)
        imshow(plane_mask, 'InitialMagnification', 'fit');
        title('combined');
        subplot(324)
        imshow(~boundary_mask, 'InitialMagnification', 'fit');
        title('boundary');
      end

      obs_x = heightmap.X(boundary_mask);
      obs_y = heightmap.Y(boundary_mask);
      obstacle_pts = reshape([reshape(obs_x, 1, []); reshape(obs_y, 1, [])], 2, 1, []);
      theta_steps = yaw + (-pi:pi/4:pi);

      %% Get obstacles for the foot shape
      % Reorient the foot onto the terrain plane
      c = cross([0;0;1], n0);
      if norm(c) >= 1e-6
        ax = c / norm(c);
        angle = asin(norm(c));
        Rplane = axis2rotmat([ax; angle]);
      else
        Rplane = eye(3);
      end

      c_obs = iris.cspace.cspace3(obstacle_pts, collision_model.foot_shape, theta_steps, Rplane(1:2,1:2));

      %% Get obstacles for the upper body collision model
      % compute distance from the heightmap Z to the current plane
      % n' * [x;y;z] = n' * point
      % z = (n'*point - n(1:2)'*[x;y])/n(3)
      dZ = heightmap.Z - reshape((n0'*p0 - n0(1:2)'*[reshape(heightmap.X,1,[]);reshape(heightmap.Y,1,[])])/n0(3), sz);
      for j = 1:length(collision_model.body_slices.z)
        zmin = collision_model.body_slices.z(j);
        if j < length(collision_model.body_slices.z)
          zmax = collision_model.body_slices.z(j+1);
        else
          zmax = inf;
        end
        z_range_mask = dZ >= zmin & dZ <= zmax;
        boundary_mask = logical(iris.terrain_grid.component_boundary(~z_range_mask, i0));
        obs_x = heightmap.X(boundary_mask);
        obs_y = heightmap.Y(boundary_mask);
        if ~isempty(obs_x)
          obstacle_pts = reshape([reshape(obs_x, 1, []); reshape(obs_y, 1, [])], 2, 1, []);
          c_obs = cat(3, c_obs, iris.cspace.cspace3(obstacle_pts, collision_model.body_slices.xy(:,:,j), theta_steps ));
        end
      end

      %% Add a bounding polyhedron
      bounds = iris.Polyhedron.from2DVertices([heightmap.X(end,1), heightmap.X(end,end), heightmap.X(1,end), heightmap.X(1,1);
                                             heightmap.Y(end,1), heightmap.Y(end,end), heightmap.Y(1,end), heightmap.Y(1,1)]);
      % Add bounds on yaw angle
      bounds.A = [bounds.A, zeros(size(bounds.A, 1), 1);
                  zeros(2, size(bounds.A, 2)), [-1; 1]];
      bounds.b = [bounds.b; -theta_steps(1); theta_steps(end)];
      if size(options.xy_bounds.A, 2) == 2
        options.xy_bounds.A(:,end+1:3) = 0;
      end
      bounds.A = [bounds.A; options.xy_bounds.A];
      bounds.b = [bounds.b; options.xy_bounds.b];
      [A, b, C, d] = iris.inflate_region(c_obs, bounds.A, bounds.b, [x0; y0; yaw], 'require_containment', true, 'error_on_infeasible_start', options.error_on_infeasible_start);

      [A, iA] = unique(A, 'rows');
      b = b(iA);
      region = iris.TerrainRegion(A, b, C, d, p0, n0);

      if options.debug
        figure(12)
        clf
        iris.drawing.drawPolyFromVertices(iris.thirdParty.polytopes.lcon2vert(A, b)', 'r');
      end

    end

    function regions = findSafeTerrainRegions(obj, map_id, collision_model, varargin)
      p = inputParser();
      p.KeepUnmatched = true;
      p.addRequired('map_id', @isnumeric);
      p.addRequired('collision_model', @(x) isa(x, 'iris.terrain_grid.CollisionModel'));
      p.addParamValue('seeds', []);
      p.addParamValue('default_yaw', 0, @isnumeric);
      p.addParamValue('max_slope_angle', 40 * pi/180, @isnumeric);
      p.addParamValue('max_height_variation', 0.05, @isnumeric);
      p.addParamValue('xy_bounds', iris.Polyhedron(zeros(0,2),zeros(0,1)), @(x) isa(x, 'iris.Polyhedron'));
      p.addParamValue('max_num_regions', inf, @(x) x > 0);
      p.addParamValue('debug', false);
      p.parse(map_id, collision_model, varargin{:});
      options = p.Results;

      region_options = p.Unmatched;
      region_options.map_id = map_id;
      region_options.xy_bounds = options.xy_bounds;
      region_options.debug = options.debug;

      regions = iris.TerrainRegion.empty();

      foot_length = max(collision_model.foot_shape(1,:)) - min(collision_model.foot_shape(1,:));

      heightmap = obj.getHeightmap(map_id);
      slope_angles = heightmap.getSlopeAngles();
      potential_safe_grid = slope_angles < options.max_slope_angle;
      sz = size(potential_safe_grid);

      within_bounds = all(bsxfun(@le, options.xy_bounds.A * [reshape(heightmap.X, 1, []); reshape(heightmap.Y, 1, [])], options.xy_bounds.b), 1);
      potential_safe_grid = potential_safe_grid & reshape(within_bounds,sz);
      for dx = -1:1
        for dy = -1:1
          dZ = heightmap.Z(2:end-1,2:end-1) - heightmap.Z((2+dx):(end-1+dx),(2+dy):(end-1+dy));
          potential_safe_grid(2:end-1,2:end-1) = potential_safe_grid(2:end-1,2:end-1) & abs(dZ) < options.max_height_variation;
        end
      end

      excluded_grid = true(sz);
      seed_ind = 1;

      while length(regions) < options.max_num_regions

        if options.debug
          figure(3)
          clf
          imshow(potential_safe_grid, 'InitialMagnification', 'fit');
        end

        if seed_ind <= size(options.seeds, 2)
          seed = options.seeds([1,2,6], seed_ind);
          yaw = seed(3);
          seed_ind = seed_ind + 1;
          i0 = obj.xy2ind(map_id, seed(1:2));
          % seed
          % heightmap.X(i0)
          % heightmap.Y(i0)
          % disp('here')
        else
          obs_dists = iris.terrain_grid.obs_dist(potential_safe_grid);
          [max_dist, i0] = max(obs_dists(:));
          yaw = options.default_yaw;
          if max_dist < 0.5 * foot_length / mean(heightmap.resolution);
            break;
          end
        end

        try
          region = obj.getCSpaceRegionAtIndex(i0, yaw, collision_model, region_options, 'excluded_grid', excluded_grid, 'error_on_infeasible_start', true);
        catch e
          if strcmp(e.identifier, 'IRIS:InfeasibleStart')
            excluded_grid(i0) = false;
            potential_safe_grid(i0) = false;
            continue
          else
            rethrow(e);
          end
        end
        regions(end+1) = region;
        % excluded_grid(i0) = false;
        potential_safe_grid(i0) = false;

        inpoly = all(bsxfun(@minus, region.A * [reshape(heightmap.X, 1, []); reshape(heightmap.Y, 1, []); yaw + zeros(1, numel(heightmap.X))], region.b) <= foot_length / 2, 1);
        inpoly = reshape(inpoly, sz);
        potential_safe_grid = potential_safe_grid & ~inpoly;
      end
    end

    function i0 = xy2ind(obj, map_id, xy)
      heightmap = obj.getHeightmap(map_id);
      sz = size(heightmap.X);
      [~, m0] = min(abs(heightmap.Y(:,1) - xy(2)));
      [~, n0] = min(abs(heightmap.X(1,:) - xy(1)));
      i0 = sub2ind(sz, m0, n0);
    end

    function obj = addHeightmap(obj, id, heightmap)
      if isa(heightmap, 'RigidBodyHeightMapTerrain')
        [X, Y] = meshgrid(heightmap.x, heightmap.y);
        [Z, normals] = heightmap.getHeight([reshape(X, 1, []); reshape(Y, 1, [])]);
        heightmap = iris.terrain_grid.Heightmap(X, Y, reshape(Z, size(X)), normals);
      elseif isa(heightmap, 'iris.terrain_grid.Heightmap')
      else
        error('IRIS:TerrainSegmentation:BadHeightmap', 'unrecognized heightmap class: %s', class(heightmap));
      end
      obj.heightmaps(id) = heightmap;
      obj.map_id_last_added(end+1) = id;
      obj.cleanup();
    end

    function b = hasHeightmap(obj, id)
       b = obj.heightmaps.isKey(id);
     end

    function heightmap = getHeightmap(obj, id)
      heightmap = obj.heightmaps(id);
      obj.map_id_last_used(obj.map_id_last_used == id) = [];
      obj.map_id_last_used(end+1) = id;
      obj.cleanup();
    end

    function cleanup(obj)
      if length(obj.map_id_last_added) >= 2 * obj.max_stored_maps || length(obj.map_id_last_used) >= 2 * obj.max_stored_maps
        last_used = obj.map_id_last_used(max(1, end-(obj.max_stored_maps-1)):end);
        last_added = obj.map_id_last_added(max(1, end-(obj.max_stored_maps-1)):end);
        ids = num2cell(union(last_added, last_used));
        obj.heightmaps = containers.Map(ids, obj.heightmaps.values(ids));
      end
    end
  end
end

author: rdeits

repo: iris-distro

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


text:
% Test ICP scan matching algorithm

load corridor.mat;
rob = Robot;
rob.R = [-2; -1; 0];
rob.r = rob.R;
rob.u = [2; pi/1.5]; % Set robot control to move in straight line
sen = Sensor;
sen.range = 30;
sen.fov = 3*pi/2;
axis_size = 5.5;

% obs is an array of objects of type obstacle
scan_data_1 = [];
for i = 1:length(obs)
    scan_data_1_tmp = obs(i).getMeasured(sen, rob);
    scan_data_1 = [scan_data_1 scan_data_1_tmp];
end
scan_data_1 = removeDuplicateLasers(scan_data_1);
scan_data_1 = scan_data_1(2:3, :);
v = repmat(sen.noise, 1, size(scan_data_1, 2)) .* randn(2,size(scan_data_1, 2));
scan_data_1 = scan_data_1 + v;
scan_data_1 = getInvMeasurement(scan_data_1);

p1 = rob.computeTriangle('true');
p1_guess = rob.computeTriangle;

rob.R = rob.move(rob.q .* randn(2,1), 'true');
p2 = rob.computeTriangle('true');

scan_data_2 = [];
for i = 1:length(obs)
    scan_data_2_tmp = obs(i).getMeasured(sen, rob);
    scan_data_2 = [scan_data_2 scan_data_2_tmp];
end
scan_data_2 = removeDuplicateLasers(scan_data_2);
scan_data_2 = scan_data_2(2:3, :);
v = repmat(sen.noise, 1, size(scan_data_2, 2)) .* randn(2,size(scan_data_2, 2));
scan_data_2 = scan_data_2 + v;
scan_data_2 = getInvMeasurement(scan_data_2);

[R, T, corr, icp_var] = doICP(scan_data_1, scan_data_2, 1, rob.u);

da = getPiToPi(asin(R(2,1)));
rob.r(1:2) = transToGlobal(rob.r, T);
rob.r(3) = rob.r(3) + da;
p2_guess = rob.computeTriangle;

% Plot results
figure('color', 'white')
s1 = subplot(1,2,1);
h1 = plot(p1(1,:), p1(2,:), 'r-', p2(1,:), p2(2,:), 'b-', ...
    p2_guess(1,:), p2_guess(2,:), 'm--');
hold on
for i = 1:length(obs)
    obs(i).plot(s1);
end
axis square
axis([-axis_size axis_size -axis_size axis_size])
title(['Before and after wheelchair poses' char(10) 'relative to obstacles'])
legend('pose before', 'pose after', 'estimated pose after', 'obstacles')

% Rot = [cos(pi/2) -sin(pi/2); sin(pi/2) cos(pi/2)];
Rot = 1;
scan_data_2_rot = Rot*scan_data_2;
s2 = subplot(1,2,2);
h2 = plot(scan_data_1(1,:), scan_data_1(2,:), 'r+', ...
    scan_data_2_rot(1,:), scan_data_2_rot(2,:), 'b+');
axis square
axis([-axis_size axis_size -axis_size axis_size])
title(['Before and after scans' char(10) 'in local (Cartesian) frame'])
legend('scan before', 'scan after')

scan_data_2 = transToGlobal(rob.r, scan_data_2);
C = getICPCovariance(icp_var, scan_data_2(:, corr(:, 2)'));
S = 1/trace(C);
disp(['Saliency score: ' num2str(S)])

author: jaijuneja

repo: ekf-slam-matlab


text:
function isol = get_isolated(st0, id0, Mask, nt1)

% nt1 = 41;

imin = 1;
icurrent = 1;
imax = 1;
nspbatch = numel(st0);

isol = true(nspbatch, 1);

nt0 = (size(Mask,1)+1)/2;
while(icurrent<=nspbatch)
    while (st0(imax) - st0(icurrent)) <= nt1-1 && imax<nspbatch
        imax = imax+1;
    end

    while (st0(icurrent) - st0(imin)) >= nt1-1
        imin = imin+1;
    end
    for i = [imin:icurrent-1 icurrent+1:imax-1]
        if (Mask(nt0 + (st0(i) - st0(icurrent)), id0(icurrent), id0(i)))
            isol(icurrent) = false;
            break;
        end
    end

    icurrent = icurrent + 1;
end

author: cortex-lab

repo: KiloSort

save_path: github-repos/MATLAB/cortex-lab-KiloSort


text:
function [tf,te,to,sc,I,up] = tsurf_cad(F,V,varargin)
  % [tf,te,to,sc,I,up] = tsurf_cad(F,V,varargin)
  %
  % TSURF_CAD Display a triangle mesh in a CAD-like rendering style: edges,
  % shadow, floor.
  %
  % Inputs:
  %   F  #F by 3 list of face indices
  %   V  #V by 3 list of vertex positions
  % Outputs:
  %   tf  plot handle to mesh
  %   te  plot handle to sharp edges
  %   to  plot handle to boundary edges
  %   sc  plot handle to shadow
  %   I  indices into V to cut mesh along sharp edges
  %   up  callback function to reset lines based on rotation
  %

  function set_hold(v)
    if v
      hold on;
    else
      hold off;
    end
  end
  old_hold = ishold;
  %V = V*axisangle2matrix([1 0 0],-pi/2);
  N = normals(V,F);
  BC = barycenter(V,F);

  E = sharp_edges(V,F);
  E = union(sort(E,2),sort(outline(F),2),'rows');

  % cut mesh at sharp edges to get crisp normals
  [G,I] = cut_edges(F,E);
  W = V(I,:);

  if ~old_hold
    clf;
  end
  hold on;
  tf = tsurf(G,W, ...
    'FaceVertexCData',repmat(blue,size(W,1),1), ...
    'SpecularStrength',0, ...
    'DiffuseStrength',0.1, ...
    'AmbientStrength',1.0, ...
    'EdgeColor','none','FaceAlpha',0.9,fphong, ...
    varargin{:});
  te = tsurf(E(:,[1 2 2]),V,'EdgeColor',blue*0.75);
  to = tsurf([1 1 1],V,'LineWidth',2,'EdgeColor',blue*0.5);
  view(-9,34);
  axis equal;
  l = light('Position',[3 -4 5.0],'Style','infinite');
  [h,~,M,g] = add_shadow(tf,l,'Ground',[0 0 -1 min(V(:,3))-2e-3],'Fade','local','Color',[0.8 0.8 0.8]);
  % faint amient occlusion
  AO = ambient_occlusion(W,G,W,per_vertex_normals(W,G),1000);
  AO = AO*0.27;
  tf.FaceVertexCData = bsxfun(@times,tf.FaceVertexCData,1-AO);
  %set_hold(old_hold);

  % Hack so that axis doesn't change if V contains unreferenced points and
  % doesn't change.
  [BB,BF] = bounding_box(V);
  hold on;
  bf = trisurf(BF,BB(:,1),BB(:,2),BB(:,3),'CData',nan*BB(:,1),'EdgeColor','none');
  set_hold(old_hold);

  % floor board
  SV = [V ones(size(V,1),1)]*M';
  SV = bsxfun(@rdivide,SV(:,1:3),SV(:,4));

  [BB,BF] = bounding_box(SV(:,1:2));
  M = mean(bounding_box(V(:,1:2)));
  BBmM = bsxfun(@minus,BB,M);
  BB = bsxfun(@plus,sign(BBmM)*max(abs(BBmM(:))),M);
  BB = bsxfun(@plus,bsxfun(@minus,BB,mean(BB))*1.1,mean(BB));
  BB(:,3) = min(V(:,3))-4e-3;
  BB = reshape(permute(BB,[1 3 2]),[2 2 3]);
  % checkboard texture
  nx = 16;
  extent = @(BB) max(max(BB)) - min(min(BB));
  ny = ceil(extent(BB(:,:,2))/ extent(BB(:,:,1))*nx/2)*2;
  ny(isinf(ny)) = nx;
  ch = repmat(0.99-0.15*xor((mod(repmat(0:8*nx-1,8*ny,1),8*2)>7), ...
    (mod(repmat((0:8*ny-1)',1,8*nx),8*2)>7)),[1 1 3])*0.5 + 0.5;
  hold on;
  sc = surf(BB(:,:,1),BB(:,:,2),BB(:,:,3), ...
    'CData',ch,'FaceColor','texturemap', ...
    'SpecularStrength',0, 'DiffuseStrength',0, 'AmbientStrength',1);
  set_hold(old_hold);

  axis vis3d;
  camproj('persp');
  set(gca,'Visible','off');
  T = get(gca,'tightinset');
  set(gca,'position',[T(1) T(2) 1-T(1)-T(3) 1-T(2)-T(4)]);

  % Set up rotatation callbacks to hide view-dependent effects during drag
  up = @() ...
    set(to,'Faces', ...
      outline(F((sum(N.*bsxfun(@minus,BC,campos),2)<=0),:))*[1 0 0;0 1 1]) | ...
    cell2mat(cellfun(@(h) set(h,'FaceAlpha',0.5*(g*[campos 1]'<0)),h,'UniformOutput',false)) | ...
    set(sc,'FaceAlpha',1.0*(g*[campos 1]'<0));
  up();
  down = @() ...
    (exist('to','var') && ishandle(to) && isempty(set(to,'Faces',[]))) || ...
    isempty(set(rotate3d,'ActionPostCallback',[],'ActionPreCallback',[]));
  set(rotate3d,'ActionPostCallback',@(src,obj) up());
  set(rotate3d,'ActionPreCallback',@(src,obj) down());

%   for t = linspace(0,-360,60)
%    view(64+t,20);
%    up();
%    drawnow;
%    filename = 'tsurf-cad.gif';
%    figgif(filename,'nodither');
%   end
end

author: alecjacobson

repo: gptoolbox


text:
function P_matrix = Initialize_LC_P_matrix(LC_KF_config)
%Initialize_LC_P_matrix - Initializes the loosely coupled INS/GNSS KF
%error covariance matrix
%
% Software for use with "Principles of GNSS, Inertial, and Multisensor
% Integrated Navigation Systems," Second Edition.
%
% This function created 12/4/2012 by Paul Groves
%
% Inputs:
%   TC_KF_config
%     .init_att_unc           Initial attitude uncertainty per axis (rad)
%     .init_vel_unc           Initial velocity uncertainty per axis (m/s)
%     .init_pos_unc           Initial position uncertainty per axis (m)
%     .init_b_a_unc           Initial accel. bias uncertainty (m/s^2)
%     .init_b_g_unc           Initial gyro. bias uncertainty (rad/s)
%
% Outputs:
%   P_matrix              state estimation error covariance matrix

% Copyright 2012, Paul Groves
% License: BSD; see license.txt for details

% Begins

% Initialize error covariance matrix
P_matrix =  zeros(15);
P_matrix(1:3,1:3) = eye(3) * LC_KF_config.init_att_unc^2;
P_matrix(4:6,4:6) = eye(3) * LC_KF_config.init_vel_unc^2;
P_matrix(7:9,7:9) = eye(3) * LC_KF_config.init_pos_unc^2;
P_matrix(10:12,10:12) = eye(3) * LC_KF_config.init_b_a_unc^2;
P_matrix(13:15,13:15) = eye(3) * LC_KF_config.init_b_g_unc^2;

% Ends

author: awerries

repo: kalman-localization


text:
function r1mach_test ( )

%*****************************************************************************80
%
%% R1MACH_TEST reports the constants returned by R1MACH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 April 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R1MACH_TEST\n' );
  fprintf ( 1, '  R1MACH reports the value of constants associated\n' );
  fprintf ( 1, '  with real single precision computer arithmetic.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Assume that single precision numbers are stored \n' );
  fprintf ( 1, '  with a mantissa of T digits in base B, with an \n' );
  fprintf ( 1, '  exponent whose value must lie between EMIN and EMAX.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  For input arguments of 1 <= I <= 5,\n' );
  fprintf ( 1, '  R1MACH will return the following values:\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  R1MACH(1) = B^*(EMIN-1), the smallest positive magnitude.\n' );
  fprintf ( 1, '%26.16e\n', r1mach(1) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  R1MACH(2) = B^EMAX*(1-B^(-T)), the largest magnitude.\n' );
  fprintf ( 1, '%26.16e\n', r1mach(2) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  R1MACH(3) = B^(-T), the smallest relative spacing.\n' );
  fprintf ( 1, '%26.16e\n', r1mach(3) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  R1MACH(4) = B^(1-T), the largest relative spacing.\n' );
  fprintf ( 1, '%26.16e\n', r1mach(4) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  R1MACH(5) = log10(B).\n' );
  fprintf ( 1, '%26.16e\n', r1mach(5) );

  return
end


text:
%% Global Patch Collider Demo
%
% This sample trains the forest for the Global Patch Collider and stores
% output to the file "forest.yml.gz".
%
% It then finds correspondences between two images using Global Patch Collider
% and calculates error using provided ground truth flow.
%
% It will look for the file named "forest.yml.gz" with a learned forest.
% You can obtain the "forest.yml.gz" either by manually training it or by
% downloading one of the files trained on some publicly available dataset from
% here: <https://drive.google.com/open?id=0B7Hb8cfuzrIIZDFscXVYd0NBNFU>
%
% Sources:
%
% * <https://github.com/opencv/opencv_contrib/blob/3.2.0/modules/optflow/samples/gpc_train.cpp>
% * <https://github.com/opencv/opencv_contrib/blob/3.2.0/modules/optflow/samples/gpc_evaluate.cpp>
%

%% 1) Train

%%
% input training images
imgs1 = {
    fullfile(mexopencv.root(), 'test', 'RubberWhale1.png')
};
imgs2 = {
    fullfile(mexopencv.root(), 'test', 'RubberWhale2.png')
};
groundTruths = {
    fullfile(mexopencv.root(), 'test', 'RubberWhale.flo')
};
assert(isequal(numel(imgs1), numel(imgs2), numel(groundTruths)));

if exist(groundTruths{1}, 'file') ~= 2
    % attempt to download ground thruth flow from GitHub
    disp('Downloading FLO...')
    url = 'https://cdn.rawgit.com/opencv/opencv_extra/3.2.0/testdata/cv/optflow/RubberWhale.flo';
    urlwrite(url, groundTruths{1});
end

%%
% Global Patch Collider training parameters
params = {
    'MaxTreeDepth',20, ...      % Maximum tree depth to stop partitioning
    'MinNumberOfSamples',3, ... % Minimum number of samples in the node to stop partitioning
    'DescriptorType','DCT', ... % Descriptor type. Set to DCT for quality, WHT for speed
    'PrintProgress',false       % Set to false for quiet mode, set to true to print progress
};
forestDumpPath = fullfile(tempdir(), 'forest.yml.gz');

%%
% train the forest for the Global Patch Collider and save it
if exist(forestDumpPath, 'file') ~= 2
    gpc = cv.GPCForest();
    tic
    gpc.train(imgs1, imgs2, groundTruths);
    toc
    gpc.save(forestDumpPath);
end

%% 2) Evaluate

%%
% test images
fromPath = imgs1{1};
toPath = imgs2{1};
gtPath = groundTruths{1};

from = imread(fromPath);
to = imread(toPath);
flo = cv.readOpticalFlow(gtPath);

%%
% load pretrained forest
forest = cv.GPCForest();
assert(exist(forestDumpPath,'file') == 2, 'No file with a trained model');
forest.load(forestDumpPath);

%%
% find correspondences between two the images using GPC
tic
corresp = forest.findCorrespondences(from, to, 'UseOpenCL',false);
toc
fprintf('Found %d matches\n', numel(corresp));

%%
% calculate error using provided ground truth flow
gtU = flo(:,:,1);
gtV = flo(:,:,2);
a = cat(1, corresp.first);
b = cat(1, corresp.second);
ind = sub2ind(size(gtU), a(:,2), a(:,1));
gtDisplacement = [gtU(ind) gtV(ind)];
c = a + gtDisplacement;

% check for correct flow vector
mask = all(isfinite(gtDisplacement) & (gtDisplacement < 1e9), 2);
a = a(mask,:);
b = b(mask,:);
c = c(mask,:);

err = mean(sqrt(sum((b - c).^2, 2)));
fprintf('Average endpoint error = %f px.\n', err);

%%
% display flows as color images
clr = getFlowColor(b - a);
dispOut = zeros(size(from), 'single');
dispOut(:,:,3) = 1;
dispOut = cv.circle(dispOut, a, 3, 'Colors',clr, 'Thickness','Filled');
dispOut = cv.cvtColor(dispOut, 'HSV2RGB');

clr = getFlowColor(b - c, false, 32);
dispErr = zeros(size(from), 'single');
dispErr(:,:,3) = 1;
dispErr = cv.circle(dispErr, a, 3, 'Colors',clr, 'Thickness','Filled');
dispErr = cv.cvtColor(dispErr, 'HSV2RGB');

dispGT = getFlowColor([gtU(:) gtV(:)]);
dispGT = reshape(dispGT(:,1:3), [size(gtU) 3]);
dispGT = cv.cvtColor(dispGT, 'HSV2RGB');

%%
% show results
opts = {'FontScale',0.8, 'Color','k', 'LineType','AA'};
str = 'Sparse matching: Global Patch Collider';
dispOut = cv.putText(dispOut, str, [20 40], opts{:});
str = sprintf('Average EPE: %.2f', err);
dispOut = cv.putText(dispOut, str, [20 80], opts{:});
str = sprintf('Number of matches: %d', nnz(mask));
dispOut = cv.putText(dispOut, str, [20 120], opts{:});

figure(1), imshow(dispOut), title('Correspondences')
figure(2), imshow(dispErr), title('Error')
figure(3), imshow(dispGT), title('Ground Truth')

%% Helper function

function clr = getFlowColor(UV, logScale, scaleDown)
    if nargin < 2, logScale = true; end
    if nargin < 3, scaleDown = 5; end

    angle = (atan2(-UV(:,2), -UV(:,1)) + pi) * 180 / pi;
    angle(all(UV == 0, 2)) = 0;

    radius = sqrt(sum(UV.^2, 2));
    if logScale
        radius = log(radius + 1);
    end
    radius = radius ./ scaleDown;
    radius = min(radius, 1);

    clr = [angle radius];
    clr(:,3) = 1;
    clr(:,4) = 0;
end

author: kyamagu


text:
clear all;
addpath(genpath('../CoreModules'));
n_epoch=20; %training epochs
dataset_name='mnist'; %dataset name
network_name='mlp'; %network name
use_gpu=(gpuDeviceCount>0) %use gpu or not
if use_gpu
    %Requires Neural Network Toolbox to use it.
    opts.use_nntoolbox=license('test','neural_network_toolbox')
end
%function handle to prepare your data
PrepareDataFunc=@PrepareData_MNIST_MLP;
%function handle to initialize the network
NetInit=@net_init_mlp_mnist_dropout;

%automatically select learning rates
use_selective_sgd=1;
%select a new learning rate every n epochs
ssgd_search_freq=10;
learning_method=@sgd; %training method: @sgd,@rmsprop,@adagrad,@adam
opts.parameters.mom=0.9;
opts.parameters.clip=1e1;
%sgd parameter
%(unnecessary if selective-sgd is used)
sgd_lr=5e-2;

%opts.parameters.clip=1e-2;
opts.parameters.weightDecay=0;
opts.parameters.batch_size=500;

Main_Template(); %call training template

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function I = integral( f, varargin )
%INTEGRAL   Complete definite integral of DISKFUN.
%
%   I = INTEGRAL(F), returns the definite integral of a DISKFUN integrated
%   over its domain of definition.
%
%   I = INTEGRAL(F, g), returns the integral of F along the
%   curve defined by the complex-valued CHEBFUN g.
%
%   I = INTEGRAL(F, 'unitcircle') returns the integral of F along the
%   unit circle.
% See also DISKFUN/SUM2, DISKFUN/INTEGRAL2.

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

%check for empty
if isempty(f)
    I = 0;
    return
end

% Another way to do sum2(f)
if ( nargin == 1 )
    I = integral2( f );
else
    if ( ischar(varargin{1}) )
        if ( strcmpi(varargin{1},'unitcircle') )
            c = chebfun(@(t) exp(1i*t), [-pi, pi]);
        else
            error('CHEBFUN:DISKFUN:INTEGRAL:unrecognizedType',...
                'Unrecognized line integral type.  Did you mean "unitcircle"');
        end
    else
        c = varargin{1};
    end
    if ( ~isa( c, 'chebfun' ) )
        I = integral2( f, varargin{ : } );
    else                            % Line integral over a CHEBFUN
        % Make complex:
        c = c + realmin*1i;
        % Line integral:
        I = sum( feval(f, c ) .* abs( diff( c ) ) );
    end
end

end

author: chebfun


text:
function X = interaction(Xa,Xb)
% Create all interaction terms from two dummy coded design matrices.
% (See Box II in Rouder et al. 2012)
%
% INPUT
% Xa, Xb = [nrObservations nrA] and [nrObservations nrB]
%           design matrices with matching number of observation (rows)
% OUTPUT
% X = Dummy coded design matrix [nrObservations nrA*nrB]
%
nA = size(Xa,2);
nB = size(Xb,2);
Xa = repmat(Xa,[1 nB]);
Xb = repmat(Xb,[1 nA]);
ix = (repmat(1:nA:nA*nB,[nA 1]) + repmat((0:nA-1)',[1 nB]))';
Xa = Xa(:,ix(:));
X= Xa.*Xb;

end

author: fieldtrip


text:
function pstruct = tapas_condhalluc_obs2_namep(pvec)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2016 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pstruct = struct;

pstruct.be = pvec(1);
pstruct.nu = pvec(2);

return;


text:
function handles = updateStandardSpaceValue(handles, acpcCoord)
% Get the standard space popup string and set the position appropriately.
% This code needs better comments.  I did a little below, but ...
%
% Pulled out of dtiFiberUI
%
% Brian (c) VISTASOFT Team, 2012.

% The acpc coordinate was either passed in, or we get it from the position
if(~exist('acpcCoord','var'))
    acpcCoord = str2num(get(handles.editPosition, 'String'));
end

% Determine the standard space.
ssVals = get(handles.popupStandardSpace,'String');
curSs  = ssVals{get(handles.popupStandardSpace,'Value')};
curPosSs = '';
if(strcmpi(curSs,'Image'))
    % What the heck is this?
    % OK, so Bob made handles.bg a structure.  One of the fields is .mat,
    % which contains a matrix.  This matrix apparently transforms data from
    % somewhere to somewhere else.  This code, which was worse before,
    % inverts the .mat matrix from the currently selected bg.
    %
    % This should Probably be something like
    %   acpc2img = dtiGet(handles,'acpc2img')
    %   curPosSs = mrAnatomXformCoords(acpc2img);
    T = inv(handles.bg(get(handles.popupBackground,'Value')).mat);
    curPosSs = mrAnatXformCoords(T, acpcCoord);
elseif(strcmpi(curSs,'Talairach'))
    % Talairach version instead of general transform as above.
    curPosSs = mrAnatAcpc2Tal(handles.talairachScale, acpcCoord);
else
    % Check for other standard spaces
    if(isfield(handles,'t1NormParams')&&~isempty(handles.t1NormParams))
        normSs = strmatch(curSs,{handles.t1NormParams(:).name});
        if(~isempty(normSs))
            if(~isfield(handles.t1NormParams(normSs),'coordLUT') || isempty(handles.t1NormParams(normSs).coordLUT))

                lutFile = dtiGetFilenameInDT6(handles.dataFile,'lutMNI');
                if(exist(lutFile,'file') && ~isempty(lutFile))
                    disp(['Loading inverse deformation from ' lutFile '...']);
                    ni = niftiRead(lutFile);
                    handles.t1NormParams(normSs).coordLUT = squeeze(ni.data(:,:,:,1,:));
                    handles.t1NormParams(normSs).inMat = ni.qto_ijk;
                    clear ni;
                else
                    lutFile = fullfile(handles.dataDir,[curSs, '_coordLUT.nii.gz']);
                    disp('Computing inverse deformation...');
                    set(handles.editPositionTal, 'String', 'computing...');pause(0.1);
                    [defX, defY, defZ] = mrAnatInvertSn(handles.t1NormParams(normSs).sn);
                    if(max(abs(defX(:)))<127.5 && max(abs(defY(:)))<127.5 && max(abs(defZ(:)))<127.5)
                        defX(isnan(defX)) = -127; defY(isnan(defY)) = -127; defZ(isnan(defZ)) = -127;
                        handles.t1NormParams(normSs).coordLUT = int8(round(cat(4,defX,defY,defZ)));
                    else
                        defX(isnan(defX)) = -999; defY(isnan(defY)) = -999; defZ(isnan(defZ)) = -999;
                        handles.t1NormParams(normSs).coordLUT = int16(round(cat(4,defX,defY,defZ)));
                    end
                    handles.t1NormParams(normSs).inMat = inv(handles.t1NormParams(normSs).sn.VF.mat);
                    % NIFTI_INTENT_DISPVECT=1006
                    intentCode = 1006;
                    intentName = ['To' curSs];
                    % NIFTI format requires that the 4th dim is always
                    % time, so we put the deformation vector [x,y,z] in the
                    % 5th dimension.
                    tmp = reshape(handles.t1NormParams(normSs).coordLUT,[size(defX) 1 3]);
                    try
                        dtiWriteNiftiWrapper(tmp,inv(handles.t1NormParams(normSs).inMat),lutFile,1,'',intentName,intentCode);
                        dtiSetFilenameInDT6(handles.dataFile,'lutMNI',lutFile);
                    catch
                        disp('Could not save LUT transform- check permissions.');
                    end
                end
            end
            curPosSs = mrAnatXformCoords(handles.t1NormParams(normSs), acpcCoord);
        end
    end
    % Try an MNI label map
    if(isempty(curPosSs))
        % strmatch will be deprecated.  Figure out what this does and
        % replace it with strcmp
        mniSs = strmatch('MNI',{handles.t1NormParams(:).name});
        curPosSs = dtiGetBrainLabel(mrAnatXformCoords(handles.t1NormParams(mniSs), acpcCoord), curSs);
    end
end
if(isnumeric(curPosSs)), curPosSs = sprintf('%.0f, %.0f, %.0f', round(curPosSs)); end
set(handles.editPositionTal, 'String', curPosSs);

return

author: vistalab