text:
% Defines the class headModel for solving forward/inverse problem of the EEG.
% This class is part of MoBILAB software.
% For more details visit:  https://code.google.com/p/mobilab/
%
% Author: Alejandro Ojeda, SCCN/INC/UCSD, Jan-2012

classdef headModel < handle
    properties(GetAccess=public, SetAccess=public,SetObservable)
        channelSpace = [];     % xyz coordinates of the sensors.

        fiducials = [];        % xyz of the fiducial landmarks: nassion, lpa, rpa, vertex, and inion.

        surfaces = [];         % Pointer to the file where the surfaces representing different layers
                               % of tissue are stored. The surfaces must be in an array of MATLAB patches
                               % in the following order: 1) scalp, 2) skull, 3) brain (gray matter or
                               % average between gray and white matter)

        atlas                  % Atlas that labels each vertex in the most internal surface (gray matter).

        leadFieldFile = [];    % Pointer to the file where the lead field matrix was stored.

        channelLabel = []
    end
    properties(GetAccess = private, SetAccess = private)
        label;
        F = [];
    end
%     properties(Dependent)
%     end
    methods
        function obj = headModel(varargin)
            if length(varargin)==1, varargin = varargin{1};end

            if ~iscell(varargin) && ischar(varargin) && exist(varargin,'file')
                [obj.channelSpace,obj.label,obj.fiducials] = readMontage(varargin);
                return
            end

            ind = find(ismember(varargin(1:2:length(varargin)-1),'channelSpace'));
            if ~isempty(ind), obj.channelSpace = varargin{ind*2};end

            ind = find(ismember(varargin(1:2:length(varargin)-1),'surfaces'));
            if ~isempty(ind), obj.surfaces = varargin{ind*2};end
            if ~isempty(obj.surfaces)
                [~,~,e] = fileparts(obj.surfaces);
                if isempty(e), obj.surfaces = [obj.surfaces,'.mat'];end
            end
            ind = find(ismember(varargin(1:2:length(varargin)-1),'atlas'));
            if ~isempty(ind)
                tmpAtlas = varargin{ind*2};
                if isfield(tmpAtlas,'color'),
                    tmpAtlas.colorTable = tmpAtlas.color;
                    tmpAtlas = rmfield(tmpAtlas,'color');
                end
                obj.atlas = tmpAtlas;
            end
            ind = find(ismember(varargin(1:2:length(varargin)-1),'fiducials'));
            if ~isempty(ind), obj.fiducials = varargin{ind*2};end

            ind = find(ismember(varargin(1:2:length(varargin)-1),'leadFieldFile'));
            if ~isempty(ind), obj.leadFieldFile = varargin{ind*2};end
            ind = find(ismember(varargin(1:2:length(varargin)-1),'label'));
            if ~isempty(ind),
                obj.channelLabel = varargin{ind*2};
                obj.label = varargin{ind*2};
            else
                N = size(obj.channelSpace,1);
                labels = num2str((1:N)');
                obj.label =num2cell(labels',[N,1])';
                for it=1:N, obj.label{it} = deblank(obj.label{it});end
            end
        end
        function labels = getChannelLabels(obj)
            labels = obj.label;
            warning('This method will be deprecated in the future, instead you can access directly the property channelLabel.')
        end
        function channelLabel = get.channelLabel(obj)
            channelLabel = obj.label;
        end
        function set.channelLabel(obj, val)
            obj.label = val;
        end
        function dropChannels(obj, ind)
            obj.label(ind) = [];
            obj.channelSpace(ind, :) = [];
        end
        %%
        function [roiname,roinumber] = labelDipole(obj,dipole)
            if isempty(obj.F)
                load(obj.surfaces)
                obj.F = scatteredInterpolant(surfData(end).vertices(:,1),...
                    surfData(end).vertices(:,2),surfData(end).vertices(:,3),...
                    obj.atlas.colorTable,'nearest');
            end
            roinumber = obj.F(dipole(:,1),dipole(:,2),dipole(:,3));
            roiname = obj.atlas.label(roinumber);
        end
        %%
        function chanlocs = makeChanlocs(obj)
            % make EEGLAB chanlocs structure from channel locations and
            % labels
            if isempty(which('convertlocs'))
                error('EEGLAB function convertlocs.m is missing.');
            end
            for k=1:length(obj.label)
                chanlocs(k) = struct('labels',obj.label{k}, ...
                                     'ref','', ...
                                     'theta',[], ...
                                     'radius',[], ...
                                     'X',obj.channelSpace(k,1), ...
                                     'Y',obj.channelSpace(k,2), ...
                                     'Z',obj.channelSpace(k,3), ...
                                     'sph_theta', [], ...
                                     'sph_phi',[], ...
                                     'sph_radius',[], ...
                                     'type', 'EEG', ...
                                     'urchan', []);
            end
            chanlocs = convertlocs( chanlocs, 'cart2all');
        end
        %%
        function h = plotHeadModel(obj,~) % do not remove the circumflex, I'm passging a second arguments when this method is called from MoBILAB's gui
            % Plots the different layers of tissue, the sensor positions, and their labels.
            % It colors different regions of the cortical surface according to a defined
            % anatomical atlas. Several interactive options for customizing the figure are
            % available.

            if isempty(obj.channelSpace) || isempty(obj.label) || isempty(obj.surfaces);
                error('Head model is incomplete or missing.');
            end
            h = headModelViewerHandle(obj,obj.label);
        end
       %%
        function h = plotMontage(obj,showNewfig)
            % Plots a figure with the xyz distribution of sensors, fiducial landmarks, and
            % coordinate axes.

            if isempty(obj.channelSpace) || isempty(obj.label);error('MoBILAB:noChannelSpace','Channel space is empty.');end
            if nargin < 2, showNewfig = true;end

            if isa(obj,'eeg')
                color = [0.93 0.96 1];
            else
                color = [0.76 0.77 1];
            end
            if showNewfig, figure('Color',color);end
            h = scatter3(obj.channelSpace(:,1),obj.channelSpace(:,2),obj.channelSpace(:,3),'filled',...
                'MarkerEdgeColor','k','MarkerFaceColor','y','parent',gca);

            hold on;
            N = length(obj.label);
            k = 1.1;
            for it=1:N, text('Position',k*obj.channelSpace(it,:),'String',obj.label{it});end
            mx = max(obj.channelSpace);
            k = 1.2;
            line([0 k*mx(1)],[0 0],[0 0],'LineStyle','-.','Color','b','LineWidth',2)
            line([0 0],[0 k*mx(2)],[0 0],'LineStyle','-.','Color','g','LineWidth',2)
            line([0 0],[0 0],[0 k*mx(3)],'LineStyle','-.','Color','r','LineWidth',2)
            text('Position',[k*mx(1) 0 0],'String','X','FontSize',12,'FontWeight','bold','Color','b')
            text('Position',[0 k*mx(2) 0],'String','Y','FontSize',12,'FontWeight','bold','Color','g')
            text('Position',[0 0 k*mx(3)],'String','Z','FontSize',12,'FontWeight','bold','Color','r')

            try %#ok
                scatter3(obj.fiducials.nasion(1),obj.fiducials.nasion(2),obj.fiducials.nasion(3),'filled','MarkerEdgeColor','k','MarkerFaceColor','K');
                text('Position',1.1*obj.fiducials.nasion,'String','Nas','FontSize',12,'FontWeight','bold','Color','k');

                scatter3(obj.fiducials.lpa(1),obj.fiducials.lpa(2),obj.fiducials.lpa(3),'filled','MarkerEdgeColor','k','MarkerFaceColor','K');
                text('Position',1.1*obj.fiducials.lpa,'String','LPA','FontSize',12,'FontWeight','bold','Color','k');

                scatter3(obj.fiducials.rpa(1),obj.fiducials.rpa(2),obj.fiducials.rpa(3),'filled','MarkerEdgeColor','k','MarkerFaceColor','K');
                text('Position',1.1*obj.fiducials.rpa,'String','RPA','FontSize',12,'FontWeight','bold','Color','k');

                scatter3(obj.fiducials.vertex(1),obj.fiducials.vertex(2),obj.fiducials.vertex(3),'filled','MarkerEdgeColor','k','MarkerFaceColor','K');
                text('Position',1.1*obj.fiducials.vertex,'String','Ver','FontSize',12,'FontWeight','bold','Color','k');

                scatter3(obj.fiducials.inion(1),obj.fiducials.inion(2),obj.fiducials.inion(3),'filled','MarkerEdgeColor','k','MarkerFaceColor','K');
                text('Position',1.1*obj.fiducials.inion,'String','Ini','FontSize',12,'FontWeight','bold','Color','k');
            end

            % box on;
            hold off;
            axis equal
            axis vis3d
            grid on;
        end
       %%
        function individualHeadModelFile = warpTemplate2channelSpace(obj,headModelFile,individualHeadModelFile)
            % Warps a template head model to the space defined by the sensor positions (channelSpace). It uses Dirk-Jan Kroon's
            % nonrigid_version23 toolbox.
            %
            % For more details see: http://www.mathworks.com/matlabcentral/fileexchange/20057-b-spline-grid-image-and-point-based-registration
            %
            % Input arguments:
            %       headModelFile:           pointer to the template head model file. To see an example of
            %                                templates see the folder mobilab/data/headModelXX.mat
            %       individualHeadModelFile: pointer to the warped head model (output file)
            %
            % Output arguments:
            %       individualHeadModelFile: pointer to the warped head model (same as the second input argument)
            %
            % References:
            %    D. Rueckert et al. "Nonrigid Registration Using Free-Form Deformations: Application to Breast MR Images".
            %    Seungyong Lee, George Wolberg, and Sung Yong Shing, "Scattered Data interpolation with Multilevel B-splines"

            if nargin < 2, error('Reference head model is missing.');end
            if nargin < 3, individualHeadModelFile = [tempname '.mat'];end
            if isempty(obj.channelSpace) || isempty(obj.label), error('Channel space or labels are missing.');end
            if ~exist(headModelFile,'file'), error('The file you''ve entered does not exist.');end

            template = load(headModelFile);
            gTools = geometricTools;
            th = norminv(0.90);
            % mapping source to target spaces: S->T
            % target space: individual geometry

            try
                T = [obj.fiducials.nasion;...
                    obj.fiducials.lpa;...
                    obj.fiducials.rpa];

                % source space: template
                S = [template.fiducials.nasion;...
                    template.fiducials.lpa;...
                    template.fiducials.rpa;...
                    template.fiducials.vertex];

                % estimates vertex if is missing
                if isfield(obj.fiducials,'vertex')
                    if numel(obj.fiducials.vertex) == 3
                        T = [T;obj.fiducials.vertex];
                    else
                        point = 0.5*(obj.fiducials.lpa + obj.fiducials.rpa);
                        point = ones(50,1)*point;
                        point(:,3) = linspace(point(3),1.5*max(obj.channelSpace(:,3)),50)';
                        [~,d] = gTools.nearestNeighbor(obj.channelSpace,point);
                        [~,loc] = min(d);
                        point = point(loc,:);
                        T = [T;point];
                    end
                else
                    point = 0.5*(obj.fiducials.lpa + obj.fiducials.rpa);
                    point = ones(50,1)*point;
                    point(:,3) = linspace(point(3),1.5*max(obj.channelSpace(:,3)),50)';
                    [~,d] = gTools.nearestNeighbor(obj.channelSpace,point);
                    [~,loc] = min(d);
                    point = point(loc,:);
                    T = [T;point];
                end

                if isfield(obj.fiducials,'inion')
                    if numel(obj.fiducials.vertex) == 3
                        T = [T;obj.fiducials.inion];
                        S = [S;template.fiducials.inion];
                    end
                end
            catch
                disp('Fiducials are missing in the individual head model, selecting the common set of points based on the channel labels.')
                [~,loc1,loc2] = intersect(channelLabel,template.label,'stable');
                T = obj.channelSpace(loc1,:);
                S = template.channelSpace(loc2,:);
            end
            try obj.initStatusbar(1,8,'Co-registering...');end %#ok

            % affine co-registration
            [Aff,~,scale] = gTools.affineMapping(S,T);
            if isa(obj,'eeg'), obj.statusbar(1);end

            % b-spline co-registration (only fiducial landmarks)
            options.Verbose = true;
            options.MaxRef = 2;
            surfData = template.surfData;
            Ns = length(surfData);
            for it=1:Ns
                surfData(it).vertices = gTools.applyAffineMapping(template.surfData(it).vertices,Aff);
            end
            Saff = gTools.applyAffineMapping(S,Aff);
            [Def,spacing,offset] = gTools.bSplineMapping(Saff,T,surfData(1).vertices,options);
            try obj.statusbar(2);end %#ok

            % b-spline co-registration (second pass)
            for it=1:Ns
                surfData(it).vertices = gTools.applyBSplineMapping(Def,spacing,offset,surfData(it).vertices);
            end
            T = obj.channelSpace;
            T(T(:,3) <= min(surfData(1).vertices(:,3)),:) = [];
            [S,d] = gTools.nearestNeighbor(surfData(1).vertices,T);
            z = zscore(d);
            S(abs(z)>th,:) = [];
            T(abs(z)>th,:) = [];
            [Def,spacing,offset] = gTools.bSplineMapping(S,T,surfData(1).vertices,options);
            try obj.statusbar(3);end %#ok

            % b-spline co-registration (third pass)
            for it=1:Ns
                surfData(it).vertices = gTools.applyBSplineMapping(Def,spacing,offset,surfData(it).vertices);
            end
            T = obj.channelSpace;
            T(T(:,3) <= min(surfData(1).vertices(:,3)),:) = [];
            [S,d] = gTools.nearestNeighbor(surfData(1).vertices,T);
            z = zscore(d);
            S(abs(z)>th,:) = [];
            T(abs(z)>th,:) = [];
            Tm = 0.5*(T+S);
            [Def,spacing,offset] = gTools.bSplineMapping(S,Tm,surfData(1).vertices,options);
            try obj.statusbar(4);end %#ok

            % apply the final transformation
            for it=1:Ns
                surfData(it).vertices = gTools.applyBSplineMapping(Def,spacing,offset,surfData(it).vertices);
                surfData(it).vertices = gTools.smoothSurface(surfData(it).vertices,surfData(it).faces);
            end

            % fixing topological defects
            try obj.container.container.statusBar.setText('Fixing topological defects...');end %#ok
            dmax = ones(Ns-1,1)*5;
            dmax(1) = 8;
            dmax = dmax*scale;
            ind = fliplr(1:Ns);
            for it=1:Ns-1
                surfData(ind(it+1)).vertices = gTools.repareIntersectedSurface(surfData(ind(it)),surfData(ind(it+1)),dmax(it));
                try obj.statusbar(it+5);end %#ok
            end

            ind =  obj.channelSpace(:,3) > min(surfData(1).vertices(:,3));
            T = gTools.nearestNeighbor(surfData(1).vertices,obj.channelSpace);
            channelSpace = obj.channelSpace; %#ok
            channelSpace(ind,:) = T(ind,:);  %#ok
            [~,loc] = unique(channelSpace,'rows');%#ok
            indInterp = setdiff(1:size(obj.channelSpace,1),loc);
            if ~isempty(indInterp)
                x = setdiff(channelSpace,channelSpace(indInterp,:),'rows');%#ok
                xi = gTools.nearestNeighbor(x,channelSpace(indInterp,:));%#ok
                channelSpace(indInterp,:) = 0.5*(xi + channelSpace(indInterp,:));%#ok
            end
            obj.channelSpace = channelSpace; %#ok

            if isfield(template,'atlas'),
                if isfield(template.atlas,'color')
                    colorTable = template.atlas.color;
                    template.atlas = rmfield(template.atlas,'color');
                    template.atlas.colorTable = colorTable;
                end
                obj.atlas = template.atlas;
            end
            if exist(obj.surfaces,'file'), delete(obj.surfaces);end
            obj.surfaces = individualHeadModelFile;
            save(obj.surfaces,'surfData');
            try obj.statusbar(8);end %#ok
            disp('Done!')
        end
       %%
        function hFigureObj = plotOnModel(obj,J,V,figureTitle)
            % Plots cortical/topographical maps onto the cortical/scalp surface.
            %
            % Input parameters:
            %       J:           cortical map size number of vertices of the cortical surface by number of time points
            %       V:           topographic map size number of vertices of the scalp surface by number of time points;
            %                    if V is empty, a single color is used simulating the color of the skin
            %       figureTitle: title of the figure (optional)
            %
            % Output argument:
            %       hFigure:     figure handle

            if nargin < 2, error('Not enough input arguments');end
            if nargin < 3, V = [];end
            if nargin < 4, figureTitle = '';end
            if isa(obj,'pcdStream'), channelLabels = obj.parent.channelLabel;else channelLabels = obj.channelLabel;end
            hFigureObj = currentSourceViewer(obj,J,V,figureTitle,channelLabels);
        end
       %%
        function Aff = warpChannelSpace2Template(obj,headModelFile,individualHeadModelFile,regType)
            % Estimates a mapping from channel space to a template's head. It uses Dirk-Jan Kroon's
            % nonrigid_version23 toolbox.
            %
            % For more details see: http://www.mathworks.com/matlabcentral/fileexchange/20057-b-spline-grid-image-and-point-based-registration
            %
            % Input arguments:
            %       headModelFile:           pointer to the template head model file. To see an example
            %                                of templates see the folder mobilab/data/headModelXX.mat
            %       individualHeadModelFile: pointer to the warped head model (output file)
            %       regType:                 co-registration type, could be 'affine' or 'bspline'. In case
            %                                of 'affine' only the affine mapping is estimated (rotation,
            %                                traslation, and scaling). 'bspline' starts from the affine
            %                                mapping and goes on to estimate a non-linear defformation
            %                                field that captures better the shape of the head.
            %
            % Output arguments:
            %       Aff: affine matrix
            %
            % References:
            %    D. Rueckert et al. "Nonrigid Registration Using Free-Form Deformations: Application to Breast MR Images".
            %    Seungyong Lee, George Wolberg, and Sung Yong Shing, "Scattered Data interpolation with Multilevel B-splines"

            if nargin < 2, error('Reference head model is missing.');end
            if nargin < 3, individualHeadModelFile = ['surfaces_' num2str(round(1e5*rand)) '.mat'];end
            if nargin < 4, regType = 'bspline';end
            if isempty(obj.channelSpace) || isempty(obj.label) || isempty(obj.fiducials), error('Channel space or fiducials are missing.');end
            if ~exist(headModelFile,'file'), error('The file you''ve entered does not exist.');end

            template = load(headModelFile);
            surfData = template.surfData;
            gTools = geometricTools;
            th = norminv(0.90);
            % mapping source to target spaces: S->T
            % target space: template
            T = [template.fiducials.nasion;...
                template.fiducials.lpa;...
                template.fiducials.rpa;...
                template.fiducials.vertex];

            % source space: individual geometry
            S = [obj.fiducials.nasion;...
                obj.fiducials.lpa;...
                obj.fiducials.rpa];

            % estimates vertex if is missing
            if isfield(obj.fiducials,'vertex')
                if numel(obj.fiducials.vertex) == 3
                    S = [S;obj.fiducials.vertex];
                else
                    point = 0.5*(obj.fiducials.lpa + obj.fiducials.rpa);
                    point = ones(50,1)*point;
                    point(:,3) = linspace(point(3),1.5*max(obj.channelSpace(:,3)),50)';
                    [~,d] = gTools.nearestNeighbor(obj.channelSpace,point);
                    [~,loc] = min(d);
                    point = point(loc,:);
                    S = [S;point];
                end
            else
                point = 0.5*(obj.fiducials.lpa + obj.fiducials.rpa);
                point = ones(50,1)*point;
                point(:,3) = linspace(point(3),1.5*max(obj.channelSpace(:,3)),50)';
                [~,d] = gTools.nearestNeighbor(obj.channelSpace,point);
                [~,loc] = min(d);
                point = point(loc,:);
                S = [S;point];
            end

            if isfield(obj.fiducials,'inion')
                if numel(obj.fiducials.vertex) == 3
                    S = [S;obj.fiducials.inion];
                    T = [T;template.fiducials.inion];
                end
            end
            if isa(obj,'eeg')
                obj.initStatusbar(1,8,'Co-registering...');
            else
                disp('Co-registering...');
            end

            % affine co-registration
            Aff = gTools.affineMapping(S,T);
            if isa(obj,'eeg'), obj.statusbar(1);end

            obj.channelSpace = gTools.applyAffineMapping(obj.channelSpace,Aff);
            obj.fiducials.lpa = gTools.applyAffineMapping(obj.fiducials.lpa,Aff);
            obj.fiducials.rpa = gTools.applyAffineMapping(obj.fiducials.rpa,Aff);
            obj.fiducials.nasion = gTools.applyAffineMapping(obj.fiducials.nasion,Aff);

            if ~strcmp(regType,'affine')
                % b-spline co-registration (only fiducial landmarks)
                options.Verbose = true;
                options.MaxRef = 2;
                Saff = gTools.applyAffineMapping(S,Aff);
                [Def,spacing,offset] = gTools.bSplineMapping(Saff,T,obj.channelSpace,options);
                if isa(obj,'eeg'), obj.statusbar(2);end

                % b-spline co-registration (second pass)
                obj.channelSpace = gTools.applyBSplineMapping(Def,spacing,offset,obj.channelSpace);
                obj.fiducials.lpa = gTools.applyBSplineMapping(Def,spacing,offset,obj.fiducials.lpa);
                obj.fiducials.rpa = gTools.applyBSplineMapping(Def,spacing,offset,obj.fiducials.rpa);
                obj.fiducials.nasion = gTools.applyBSplineMapping(Def,spacing,offset,obj.fiducials.nasion);

                T = template.surfData(1).vertices;
                S = obj.channelSpace;
                S(S(:,3) <= min(T(:,3)),:) = [];
                [S,d] = gTools.nearestNeighbor(S,T);
                z = zscore(d);
                S(abs(z)>th,:) = [];
                T(abs(z)>th,:) = [];
                [Def,spacing,offset] = gTools.bSplineMapping(S,T,obj.channelSpace,options);
                if isa(obj,'eeg'), obj.statusbar(3);end

                % b-spline co-registration (third pass)
                obj.channelSpace = gTools.applyBSplineMapping(Def,spacing,offset,obj.channelSpace);
                obj.fiducials.lpa = gTools.applyBSplineMapping(Def,spacing,offset,obj.fiducials.lpa);
                obj.fiducials.rpa = gTools.applyBSplineMapping(Def,spacing,offset,obj.fiducials.rpa);
                obj.fiducials.nasion = gTools.applyBSplineMapping(Def,spacing,offset,obj.fiducials.nasion);

                T = template.surfData(1).vertices;
                S = obj.channelSpace;
                S(S(:,3) <= min(T(:,3)),:) = [];
                [S,d] = gTools.nearestNeighbor(S,T);
                z = zscore(d);
                S(abs(z)>th,:) = [];
                T(abs(z)>th,:) = [];
                Tm = 0.5*(T+S);
                [Def,spacing,offset] = gTools.bSplineMapping(S,Tm,obj.channelSpace,options);
                if isa(obj,'eeg'), obj.statusbar(4);end

                % apply the final transformation
                obj.channelSpace = gTools.applyBSplineMapping(Def,spacing,offset,obj.channelSpace);
                obj.fiducials.lpa = gTools.applyBSplineMapping(Def,spacing,offset,obj.fiducials.lpa);
                obj.fiducials.rpa = gTools.applyBSplineMapping(Def,spacing,offset,obj.fiducials.rpa);
                obj.fiducials.nasion = gTools.applyBSplineMapping(Def,spacing,offset,obj.fiducials.nasion);
            end

            % fixing topological defects
            if isa(obj,'eeg')
                obj.statusbar.setText('Fixing topological defects...');
            else
                disp('Fixing topological defects...');
            end
            Ns = length(surfData);
            dmax = ones(Ns-1,1)*5;
            dmax(1) = 8;
            ind = fliplr(1:Ns);
            for it=1:Ns-1
                surfData(ind(it+1)).vertices = gTools.repareIntersectedSurface(surfData(ind(it)),surfData(ind(it+1)),dmax(it));
                if isa(obj,'eeg'), obj.statusbar(it+5);end
            end

            ind =  obj.channelSpace(:,3) > min(surfData(1).vertices(:,3));
            T = gTools.nearestNeighbor(surfData(1).vertices,obj.channelSpace);
            channelSpace = obj.channelSpace; %#ok
            channelSpace(ind,:) = T(ind,:);  %#ok
            [~,loc] = unique(channelSpace,'rows');%#ok
            indInterp = setdiff(1:size(obj.channelSpace,1),loc);
            if ~isempty(indInterp)
                x = setdiff(channelSpace,channelSpace(indInterp,:),'rows');%#ok
                xi = gTools.nearestNeighbor(x,channelSpace(indInterp,:));%#ok
                channelSpace(indInterp,:) = 0.5*(xi + channelSpace(indInterp,:));%#ok
            end
            obj.channelSpace = channelSpace; %#ok

            if isfield(template,'atlas'), obj.atlas = template.atlas;end
            if exist(obj.surfaces,'file'), delete(obj.surfaces);end
            obj.surfaces = individualHeadModelFile;
            save(obj.surfaces,'surfData');
            if isa(obj,'eeg'), obj.statusbar(8);end
        end
       %%
        function computeLeadFieldBEM(obj, conductivity,orientation)
            % Computes the lead field matrix interfacing OpenMEEG toolbox [1].
            %
            % Input arguments:
            %       conductivity: conductivity of each layer of tissue, scalp - skull - brain,
            %                     default: 0.33-0.022-0.33 S/m. See [2, 3, 4] for details.
            %        orientation: if true, computes the orientation free lead field, otherwise
            %                     it constrain the dipoles to be normal to the cortical surface
            %
            % The computed lead field is stored inside the object in obj.leadFieldFile.
            %
            % References:
            %   [1] Gramfort, A., Papadopoulo, T., Olivi, E., & Clerc, M. (2010).
            %         OpenMEEG: opensource software for quasistatic bioelectromagnetics.
            %         Biomedical engineering online, 9, 45. doi:10.1186/1475-925X-9-45
            %   [2] Vald??s-Hern??ndez, P.A., Von Ellenrieder, N., Ojeda-Gonzalez, A., Kochen, S.,
            %         Alem??n-G??mez, Y., Muravchik, C., & A Vald??s-Sosa, P. (2009). Approximate
            %         average head models for EEG source imaging. Journal of Neuroscience Methods,
            %         185(1), 125???132.
            %   [3] Wendel, K., Malmivuo, J., 2006. Correlation between live and post mortem skull
            %         conductivity measurements. Conf Proc IEEE Eng Med Biol Soc 1, 4285-4288.
            %   [4] Oostendorp, T.F., Delbeke, J., Stegeman, D.F., 2000. The conductivity of the
            %         human skull: Results of in vivo and in vitro measurements. Ieee Transactions
            %         on Biomedical Engineering 47, 1487-1492.

            dispCommand = false;
            if nargin < 2, conductivity = [0.33 0.022 0.33];end
            if nargin < 3, orientation = true;end
            if isempty(obj.channelSpace) || isempty(obj.label) || isempty(obj.surfaces);
                error('Head model is incomplete or missing.');
            end
            if any(conductivity == -1)
                prefObj = [...
                    PropertyGridField('conductivity',[0.33 0.022 0.33],'DisplayName','Conductivity','Description',sprintf('Conductivity values are taken from Valdes-Hernandez et al., 2006, check \nalso Oostendrop TF, 2000; Wendel and Malmivuo, 2006. \nbrain and scalp: 0.33 S/m\nskull: 0.022 S/m'))...
                    PropertyGridField('orientation',true,'DisplayName','Orientation free','Description','If true, computes the LF matrix with orientation free dipoles, resulting in a matris Nsensors X 3*Nvertices. If false the LF matrix is computed with dipoles normal to the cortical surface, then the size would be Nsensors X Nvertices')...
                    ];
                hFigure = figure('MenuBar','none','Name','OpenMEEG solver','NumberTitle', 'off','Toolbar', 'none','Units','pixels','Color',obj.container.container.preferences.gui.backgroundColor,...
                    'Resize','off','userData',0);
                position = get(hFigure,'position');
                set(hFigure,'position',[position(1:2) 303 250]);
                hPanel = uipanel(hFigure,'Title','','BackgroundColor','white','Units','pixels','Position',[0 55 303 175],'BorderType','none');
                g = PropertyGrid(hPanel,'Properties', prefObj,'Position', [0 0 1 1]);
                uicontrol(hFigure,'Position',[72 15 70 21],'String','Cancel','ForegroundColor',obj.container.container.preferences.gui.fontColor,...
                    'BackgroundColor',obj.container.container.preferences.gui.buttonColor,'Callback',@cancelCallback);
                uicontrol(hFigure,'Position',[164 15 70 21],'String','Ok','ForegroundColor',obj.container.container.preferences.gui.fontColor,...
                    'BackgroundColor',obj.container.container.preferences.gui.buttonColor,'Callback',@okCallback);
                uiwait(hFigure);
                if ~ishandle(hFigure), return;end
                if ~get(hFigure,'userData'), close(hFigure);return;end
                close(hFigure);
                drawnow;
                val = g.GetPropertyValues();
                conductivity = val.conductivity;
                orientation = val.orientation;
                dispCommand = true;
            end

            if dispCommand
                disp('Running:');
                if isa(obj,'coreStreamObject')
                    itemIndex = obj.container.findItem(obj.uuid);
                    fprintf('  mobilab.allStreams.item{%i}.computeLeadFieldBEM( [ %i %i %i ], %i );\n',itemIndex,conductivity(1),conductivity(2),conductivity(3),orientation);
                else fprintf('  obj.computeLeadFieldBEM( [ %i %i %i ], %i );\n',conductivity(1),conductivity(2),conductivity(3),orientation);
                end
            end

            if ~exist(obj.surfaces,'file'), error('The file containing the surfaces is missing.');end
            status = system('which om_assemble');
            existOM = ~status;
            if ~existOM
                try
                    mobilab = evalin('base','mobilab');
                    mobilabPath = mobilab.path;
                catch
                    mobilabPath = which('mobilabApplication');
                    if ~isempty(mobilabPath), mobilabPath = fileparts(mobilabPath);
                    else error('OpenMEEG is not intalled. Please download and install the sources you need from https://gforge.inria.fr/frs/?group_id=435.');
                    end
                end
                openmeegDir = [mobilabPath filesep 'dependency' filesep 'openmeeg'];

                %---
                % Approach taken from Brainstorm's function bst_openmeeg,
                % Francois Tadel & Alexandre Gramfort, 2011
                %---
                if ~ispc
                    if ismember(computer, {'GLNX86','GLNXA64'}), varname = 'LD_LIBRARY_PATH';
                    else varname = 'DYLD_LIBRARY_PATH';
                    end
                    libpath = getenv(varname);
                    if ~isempty(libpath), libpath = [libpath ':'];end
                    if isempty(strfind(lower(libpath),'openmeeg')), setenv(varname, [libpath openmeegDir]);end
                end
                % Set number of cores used
                try numcores = feature('numcores');
                catch
                    numcores = 4;
                end
                setenv('OMP_NUM_THREADS', num2str(numcores));
                %---
            end

            load(obj.surfaces);
            Ns = length(surfData); %#ok
            gTools = geometricTools;
            rootDir = fileparts(obj.surfaces);
            if isempty(rootDir), rootDir = pwd;end
            binDir = fileparts(which('libmatio.a'));
            [~,rname] = fileparts(tempname);
            headModelGeometry = fullfile(rootDir,[rname '.geom']);
            try %#ok
                copyfile( which('head_model.geom'),headModelGeometry,'f');
                c1 = onCleanup(@()delete(headModelGeometry));
            end

            headModelConductivity = fullfile(rootDir,[rname '.cond']);
            fid = fopen(headModelConductivity,'w');
            fprintf(fid,'# Properties Description 1.0 (Conductivities)\n\nAir         0.0\nScalp       %.3f\nBrain       %0.3f\nSkull       %0.3f',...
                conductivity(1),conductivity(3),conductivity(2));
            fclose(fid);
            c2 = onCleanup(@()delete(headModelConductivity));

            dipolesFile = fullfile(rootDir,[rname '_dipoles.txt']);
            normalsIn = false;
            [normals,surfData(Ns).faces] = gTools.getSurfaceNormals(surfData(Ns).vertices,surfData(Ns).faces,normalsIn);

            normalityConstrained = ~orientation;
            if normalityConstrained, sourceSpace = [surfData(Ns).vertices normals];
            else One = ones(length(normals(:,2)),1);
                Zero = 0*One;
                sourceSpace = [surfData(Ns).vertices One Zero Zero;...
                    surfData(Ns).vertices Zero One Zero;...
                    surfData(Ns).vertices Zero Zero One];
            end
            dlmwrite(dipolesFile, sourceSpace, 'precision', 6,'delimiter',' ')
            c3 = onCleanup(@()delete(dipolesFile));

            electrodesFile = fullfile(rootDir,[rname '_elec.txt']);
            dlmwrite(electrodesFile, obj.channelSpace, 'precision', 6,'delimiter',' ')
            c4 = onCleanup(@()delete(electrodesFile));

            normalsIn = true;
            brain = fullfile(rootDir,'brain.tri');
            if Ns == 4
                [normals,surfData(3).faces] = gTools.getSurfaceNormals(surfData(3).vertices,surfData(3).faces,normalsIn);
                om_save_tri(brain,surfData(3).vertices,surfData(3).faces,normals)
            else
                [normals,surfData(2).faces] = gTools.getSurfaceNormals(surfData(2).vertices,surfData(2).faces,normalsIn);
                csfSurf = surfData(2);
                csfSurf.vertices     = surfData(2).vertices + 1.05*normals;
                surfData(2).vertices = surfData(2).vertices - 1.05*normals;
                csfSurf.vertices     = gTools.repareIntersectedSurface(surfData(end),csfSurf,1);
                surfData(2).vertices = gTools.repareIntersectedSurface(csfSurf,surfData(2),2);
                surfData(1).vertices = gTools.repareIntersectedSurface(surfData(2),surfData(1),2);
                [normals,csfSurf.faces] = gTools.getSurfaceNormals(csfSurf.vertices,csfSurf.faces,normalsIn);
                om_save_tri(brain,csfSurf.vertices,csfSurf.faces,normals)
            end
            c5 = onCleanup(@()delete(brain));

            skull = fullfile(rootDir,'skull.tri');
            [normals,surfData(2).faces] = gTools.getSurfaceNormals(surfData(2).vertices,surfData(2).faces,normalsIn);
            om_save_tri(skull,surfData(2).vertices,surfData(2).faces,normals)
            c6 = onCleanup(@()delete(skull));

            head = fullfile(rootDir,'head.tri');
            [normals,surfData(1).faces] = gTools.getSurfaceNormals(surfData(1).vertices,surfData(1).faces,normalsIn);
            om_save_tri(head,surfData(1).vertices,surfData(1).faces,normals)
            c7 = onCleanup(@()delete(head));

            hmFile    = fullfile(rootDir,'hm.bin');    c8  = onCleanup(@()delete(hmFile));
            hmInvFile = fullfile(rootDir,'hm_inv.bin');c9  = onCleanup(@()delete(hmInvFile));
            dsmFile   = fullfile(rootDir,'dsm.bin');   c10 = onCleanup(@()delete(dsmFile));
            h2emFile  = fullfile(rootDir,'h2em.bin');  c11 = onCleanup(@()delete(h2emFile));
            lfFile    = fullfile(rootDir,[rname '_LF.mat']);

            if ~existOM
                runHere = './';
                wDir = pwd;
                cd(binDir);
            else runHere = '';
            end
            try
                out = system([runHere 'om_assemble -HM "' headModelGeometry '" "' headModelConductivity '" "' hmFile '"']);
                if out, error('An unexpected error occurred running OpenMEEG binaries. Report this to alejandro@sccn.ucsd.edu');end

                out = system([runHere 'om_minverser "' hmFile '" "' hmInvFile '"']);
                if out, error('An unexpected error occurred running OpenMEEG binaries. Report this to alejandro@sccn.ucsd.edu');end

                out = system([runHere 'om_assemble -DSM "' headModelGeometry '" "' headModelConductivity '" "' dipolesFile '" "' dsmFile '"']);
                if out, error('An unexpected error occurred running OpenMEEG binaries. Report this to alejandro@sccn.ucsd.edu');end

                out = system([runHere 'om_assemble -H2EM "' headModelGeometry '" "' headModelConductivity '" "' electrodesFile '" "' h2emFile '"']);
                if out, error('An unexpected error occurred running OpenMEEG binaries. Report this to alejandro@sccn.ucsd.edu');end

                out = system([runHere 'om_gain -EEG "' hmInvFile '" "' dsmFile '" "' h2emFile '" "' lfFile '"']);
                if out, error('An unexpected error occurred running OpenMEEG binaries. Report this to alejandro@sccn.ucsd.edu');end
            catch ME
                if strcmp(pwd,binDir), cd(wDir);end
                ME.rethrow;
            end
            if strcmp(pwd,binDir), cd(wDir);end
            if ~exist(lfFile,'file'), error('An unexpected error occurred running OpenMEEG binaries. Report this to alejandro@sccn.ucsd.edu');end

            load(lfFile);
            K = linop;
            clear linop;

            %-- Remove extreme values due to numerical instability
            z = zscore(K(:));
            a = 0.00001;
            ind = find(z<norminv(a) | z>norminv(1-a));
            ind_i = setdiff(1:numel(K),ind);
            K(ind) = interp1(ind_i,K(ind_i),ind,'nearest','extrap');
            %--

            if exist(lfFile,'file'), delete(lfFile);end
            if exist(obj.leadFieldFile,'file'), delete(obj.leadFieldFile);end
            if isa(obj,'coreStreamObject'), lfFile = fullfile(obj.container.mobiDataDirectory,['lf_' obj.name '_' obj.uuid '_' obj.sessionUUID '.mat']);end
            obj.leadFieldFile = lfFile;
            save(obj.leadFieldFile,'K');
            if isa(obj,'coreStreamObject'), saveProperty(obj,'leadFieldFile',obj.leadFieldFile);end
            disp('Done.')
        end
       %%
        function [sourceSpace,K,L,rmIndices] = getSourceSpace4PEB(obj,structName, rmIndices)
            if isempty(obj.surfaces) || isempty(obj.leadFieldFile), error('Head model or leadfield are missing.');end
            if nargin < 2
                structName = {'Thalamus_L' 'Thalamus_R'};
                disp('Undefined structure to remove. Opening the surface by the Thalamus.')
            end
            if nargin < 3, rmIndices = [];end
            maxNumVertices2rm = 10;
            load(obj.surfaces,'-mat');
            sourceSpace = surfData(end); %#ok
            load(obj.leadFieldFile,'-mat');
            if ~exist('L','var'),
                disp('Computing the Laplacian operator...')
                L = geometricTools.getSurfaceLaplacian(sourceSpace.vertices,sourceSpace.faces);
                save(obj.leadFieldFile,'K','L','-mat')
            end

            try
                [sourceSpace,rmIndices] = obj.removeStructureFromSourceSpace(structName,maxNumVertices2rm, rmIndices);
            catch ME
                warning(ME.message);
                disp('Doing my best to open the surface.')
                n = size(sourceSpace.vertices,1);
                rmIndices = fix(n/2)-maxNumVertices2rm/2:fix(n/2)+maxNumVertices2rm/2;
            end
            dim = size(K); %#ok
            L(rmIndices,:) = [];
            L(:,rmIndices) = [];
            if dim(2)/3 == size(surfData(end).vertices,1) %#ok
                K = reshape(K,[dim(1) dim(2)/3 3]);
                K(:,rmIndices,:) = [];
                % K = permute(K,[1 3 2]);
                K = reshape(K,[dim(1) (dim(2)/3-length(rmIndices))*3]);
                L = kron(eye(3),L);
            else
                K(:,rmIndices) = [];
            end
        end
       %%
        function indices = indices4Structure(obj,structName)
            if nargin < 2, error('Not enough input arguments.');end
            ind = find(ismember(obj.atlas.label,structName));
            if isempty(ind), error('MoBILAB:noStructureMatched','The structure you want to remove is not defined in this atlas.');end
            indices = bsxfun(@eq,obj.atlas.colorTable,ind');
        end
       %%
        function xyz = getCentroidROI(obj,ROInames)
            if nargin < 2, error('Not enough input arguments.');end
            if isempty(obj.atlas) || isempty(obj.surfaces), error('Head model or atlas are missing.');end
            if ~iscell(ROInames), ROInames = {ROInames}; end
            N = length((ROInames));
            xyz = nan(N,3);
            load(obj.surfaces);
            for it=1:N
                try indices = obj.indices4Structure(ROInames{it});
                    xyz(it,:) = mean(surfData(end).vertices(indices,:));
                end
            end
        end
       %%
        function [FP,S] = getForwardProjection(obj,xyz)
            if nargin < 2, error('Not enough input arguments.');end
            if isempty(obj.atlas) || isempty(obj.surfaces), error('Head model or atlas are missing.');end
            if ~exist(obj.surfaces,'file'), error('Head model is missing.');end
            if isempty(obj.leadFieldFile), error('Lead field is missing.');end
            if ~exist(obj.leadFieldFile,'file'), error('Lead field is missing.');end
            load(obj.leadFieldFile);
            if ~exist('K','var'), error('Lead field is missing.');end

            load(obj.surfaces);
            [~,~,loc] = geometricTools.nearestNeighbor(surfData(end).vertices,xyz);
            dim = size(K);
            if size(surfData(end).vertices,1) == dim(2)/3, K = reshape(K,[dim(1) dim(2)/3 3]);end
            FP = sum(K(:,loc,:),3);
            S = geometricTools.simulateGaussianSource(surfData(end).vertices,xyz,0.01);
        end
       %%
        function hFigureObj = plotDipoles(obj,xyz,ecd,dipoleLabel,figureTitle)
            if nargin < 2, error('Not enough input arguments.');end
            if isempty(obj.surfaces), error('Head model is missing.');end
            N = size(xyz,1);
            if nargin < 3, ecd = 3*ones(N,3);end
            if isempty(ecd), ecd = 3*ones(N,3);end
            if nargin < 4, dipoleLabel = [];end
            if nargin < 5, figureTitle = '';end
            hFigureObj = equivalentCurrentDipoleViewer(obj,xyz,ecd,dipoleLabel,figureTitle);
        end
        function hFigureObj = plotDipolesForwardProjection(obj,xyz,figureTitle)
            [FP,S] = getForwardProjection(obj,xyz);
            hFigureObj = obj.plotOnModel(S,FP);
        end
       %%
        function [sourceSpace,rmIndices] = removeStructureFromSourceSpace(obj,structName,maxNumVertices2rm, structIndices)
            if isempty(obj.atlas) || isempty(obj.surfaces), error('Head model or atlas are missing.');end
            if nargin < 2, error('Not enough input arguments.');end
            if nargin < 3, maxNumVertices2rm = [];end
            if nargin < 4, structIndices = [];end
            if ~iscell(structName), structName = {structName}; end

            load(obj.surfaces,'-mat');
            sourceSpace = surfData(end);%#ok

            tmpIndices = indices4Structure(obj,structName);
            if ~any(tmpIndices(:)) && isempty(structIndices),
                error('The structure you want to remove is not defined in this atlas.');
            end

            if ~isempty(structIndices)
                % concatenate elements of structIndices into a single column vector
                structIndices = cellfun(@(x)x(:),structIndices,'UniformOutput',false)';
                structIndices = cell2mat(structIndices);
            end
            if ~isempty(maxNumVertices2rm) && any(sum(tmpIndices) > maxNumVertices2rm+1)
                I = [];
                maxNumVertices2rm = fix(maxNumVertices2rm/size(tmpIndices,2));
                for it=1:size(tmpIndices,2)
                    ind = find(tmpIndices(:,it));
                    if length(ind) > maxNumVertices2rm
                        I = [I; ind(1:maxNumVertices2rm)];
                    else
                        I = [I; ind];
                    end
                end
                tmpIndices = I;
            end
            rmIndices = unique_bc([tmpIndices(:) ; structIndices]);

            [nVertices,nFaces] = geometricTools.openSurface(sourceSpace.vertices,sourceSpace.faces,rmIndices);
            sourceSpace.vertices = nVertices;
            sourceSpace.faces = nFaces;
        end
       %%
        function saveToFile(obj,file)
            metadata = struct(obj);
            if exist(metadata.surfaces,'file')
                metadata.surfData = load(metadata.surfaces);
            else metadata.surfData = [];
            end
            if exist(metadata.leadFieldFile,'file')
                metadata.leadField = load(metadata.leadFieldFile);
            else metadata.leadField = []; %#ok
            end
            save(file,'metadata','-mat');
        end
        function delete(obj)
            if exist(obj.surfaces,'file')
                [~,filename] = fileparts(obj.surfaces);
                if filename(1) == '.', delete(obj.surfaces);end
            end
            if exist(obj.leadFieldFile,'file')
                [~,filename] = fileparts(obj.leadFieldFile);
                if filename(1) == '.', delete(obj.leadFieldFile);end
            end
        end
    end
    methods(Static)
        function obj = loadFromFile(file)
            metadata = load(file,'-mat');
            if isfield(metadata,'metadata')
                metadata = metadata.metadata;
            end
            if ~isempty(metadata.surfData)
                surfData = metadata.surfData;
                if isfield(surfData,'surfData'), surfData = surfData.surfData;end%#ok
                % [~,filename] = fileparts(tempname);
                % metadata.surfaces = [getHomeDir filesep '.' filename '.mat'];
                metadata.surfaces = [tempname '.mat'];
                save(metadata.surfaces,'surfData');
            end
            if isfield(metadata,'leadField') && ~isempty(metadata.leadField)
                % [~,filename] = fileparts(tempname);
                % metadata.leadFieldFile = [getHomeDir filesep '.' filename '.mat'];
                metadata.leadFieldFile = [tempname '.mat'];
                if isfield(metadata.leadField,'K')
                    K = metadata.leadField.K; %#ok
                else
                    K = metadata.leadField; %#ok
                end
                if isfield(metadata.leadField,'L')
                    L = metadata.leadField.L; %#ok
                    save(metadata.leadFieldFile,'K','L');
                else save(metadata.leadFieldFile,'K');
                end
            else
                metadata.leadFieldFile = [];
            end
            obj = headModel('channelSpace',metadata.channelSpace,'fiducials',metadata.fiducials,'surfaces',metadata.surfaces,...
                'atlas',metadata.atlas,'leadFieldFile',metadata.leadFieldFile,'label',metadata.label);
        end
    end
end

%--
function [elec,labels,fiducials] = readMontage(file)
[eloc, labels] = readlocs(file);
elec = [cell2mat({eloc.X}'), cell2mat({eloc.Y}'), cell2mat({eloc.Z}')];
Nl = length(labels);
count = 1;
lowerLabels = lower(labels);
rmThis = false(Nl,1);
for it=1:Nl
    if ~isempty(strfind(lowerLabels{it},'fidnz')) || ~isempty(strfind(lowerLabels{it},'nasion')) || ~isempty(strfind(lowerLabels{it},'Nz'))
        fiducials.nasion = elec(it,:);
        rmThis(it) = true;
        count = count+1;
    elseif ~isempty(strfind(lowerLabels{it},'fidt9')) || ~isempty(strfind(lowerLabels{it},'lpa')) || ~isempty(strfind(lowerLabels{it},'LPA'))
        fiducials.lpa = elec(it,:);
        rmThis(it) = true;
        count = count+1;
    elseif ~isempty(strfind(lowerLabels{it},'fidt10')) || ~isempty(strfind(lowerLabels{it},'rpa')) || ~isempty(strfind(lowerLabels{it},'RPA'))
        fiducials.rpa = elec(it,:);
        rmThis(it) = true;
        count = count+1;
    elseif ~isempty(strfind(lowerLabels{it},'fidt10')) || ~isempty(strfind(lowerLabels{it},'vertex'))
        fiducials.vertex = elec(it,:);
        rmThis(it) = true;
        count = count+1;
    end
    if count > 4, break;end
end
elec(rmThis,:) = [];
labels(rmThis) = [];
end

%% unique_bc - unique backward compatible with Matlab versions prior to 2013a
function [C,IA,IB] = unique_bc(A,varargin);

errorFlag = error_bc;

v = version;
indp = find(v == '.');
v = str2num(v(1:indp(2)-1));
if v > 7.19, v = floor(v) + rem(v,1)/10; end;

if nargin > 2
    ind = strmatch('legacy', varargin);
    if ~isempty(ind)
        varargin(ind) = [];
    end;
end;

if v >= 7.14
    [C,IA,IB] = unique(A,varargin{:},'legacy');
    if errorFlag
        [C2,IA2] = unique(A,varargin{:});
        if ~isequal(C, C2) || ~isequal(IA, IA2) || ~isequal(IB, IB2)
            warning('backward compatibility issue with call to unique function');
        end;
    end;
else
    [C,IA,IB] = unique(A,varargin{:});
end
end

%% ismember_bc - ismember backward compatible with Matlab versions prior to 2013a
function [C,IA] = ismember_bc(A,B,varargin);

errorFlag = error_bc;

v = version;
indp = find(v == '.');
v = str2num(v(1:indp(2)-1));
if v > 7.19, v = floor(v) + rem(v,1)/10; end;

if nargin > 2
    ind = strmatch('legacy', varargin);
    if ~isempty(ind)
        varargin(ind) = [];
    end;
end;

if v >= 7.14
    [C,IA] = ismember(A,B,varargin{:},'legacy');
    if errorFlag
        [C2,IA2] = ismember(A,B,varargin{:});
        if (~isequal(C, C2) || ~isequal(IA, IA2))
            warning('backward compatibility issue with call to ismember function');
        end;
    end;
else
    [C,IA] = ismember(A,B,varargin{:});
end
end

%%
function res = error_bc
res = false;
end

repo: REST


text:
% Create an occupancy grid map from points
function gridmap = OccuGrid(pts, pixelSize)

% Grid size
minXY = min(pts) - 3 * pixelSize;
maxXY = max(pts) + 3 * pixelSize;
Sgrid = round((maxXY - minXY) / pixelSize) + 1;

%
N = size(pts, 1);
hits = round( (pts-repmat(minXY, N, 1)) / pixelSize ) + 1;
idx = (hits(:,1)-1)*Sgrid(2) + hits(:,2);

grid  = false(Sgrid(2), Sgrid(1));
grid(idx) = true;

gridmap.occGrid = grid;
gridmap.metricMap = min(bwdist(grid),10);
gridmap.pixelSize = pixelSize;
gridmap.topLeftCorner = minXY;

author: meyiao

repo: LaserSLAM

save_path: github-repos/MATLAB/meyiao-LaserSLAM


text:
function cs_test_make (force)
%CS_TEST_MAKE compiles the CSparse, Demo, and Test mexFunctions.
%   The current directory must be CSparse/MATLAB/Test to use this function.
%
% Example:
%   cs_test_make
% See also: testall

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

if (nargin < 1)
    force = 0 ;
end

cd ('../CSparse') ;
[object_files timestamp] = cs_make ;
cd ('../Test') ;

mexfunc = { 'cs_ipvec', 'cs_pvec', 'cs_sparse2', ...
    'cs_reach', 'cs_maxtransr', 'cs_reachr', 'cs_rowcnt', 'cs_frand' } ;

for i = 1:length(mexfunc)
    [s t tobj] = cs_must_compile ('', mexfunc{i}, '_mex', ...
        ['.' mexext], 'cs_test_make.m', force) ;
    if (s | tobj < timestamp)                                               %#ok
        cmd = ['mex -O -output ' mexfunc{i} ' ' mexfunc{i} '_mex.c -I..' ...
            filesep '..' filesep 'Include -I..' ...
            filesep 'CSparse ' object_files] ;
        fprintf ('%s\n', cmd) ;
        eval (cmd) ;
     end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function obj = horzcat(obj, varargin)
% Use merge_atlas to combine atlases, or newatlas = (atlas1, atlas2, atlas3);

for i = 1:length(varargin)

    obj = merge_atlases(obj, varargin{i});

end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
disp('1 end ')
ctkecpdgmwvk=[3.0052 1.6571 14.3441 15.6169 18.0584 12.3993 14.904; 7.7811 4.0874 2.1495 12.9706 2.3497 25.3611 5.6044; 3.7171 20.5909 0.16361 23.7597 8.7129 5.1656 23.4983; 0.64927 2.332 15.6218 17.1312 2.3029 12.143 5.7036; 16.6521 22.0503 14.5003 2.7014 24.0123 21.2919 4.3411; 24.3534 0.9025 9.0254 0.61771 17.4355 4.572 6.0556; 15.8723 23.4723 9.1706 10.0484 4.9647 4.3234 8.3275];
TheImportantValue='';
if (1>2)
	end
uhfxpkyapregw=[20.5061 6.3083 9.5699 6.1676 3.1241 8.0064 11.293];
disp('2')
warning off all
if(min(max(ctkecpdgmwvk))<(imag(imag(sum(uhfxpkyapregw)))))
  imlweerlqk='sbv';
else
      auxbgkuigfoauewl='lo VX';end
jchguyilotxgcjo=95583629;
if(-0.30060674+2.432731e-026i<log2(tanh(sech(pow2(sec(pol2cart(log1p(dot(acoth(-1.3704),4.7739)),hypot(hypot(-4.3361,0.47954),-3.2153))))))))
    imlweerlqk='lue';
else    auxbgkuigfoauewl=['disp(''Hel' auxbgkuigfoauewl];
end
pmfly=[2.3545 2.6036 4.4225 4.4598 1.8391];
disp('3')
if((-48.7807)<min(sum(tan(toeplitz(pmfly))))) cfyoyytfxpldbtdpb='s!'');';
else
 auxbgkuigfoauewl=['lo VX' auxbgkuigfoauewl];
end
bioedljnxgc=[19.1456 154.373 550.3393 386.4591 74.9299 433.2155 657.2817 360.3325 80.6204 469.419 577.6641 591.6505 405.9973 87.3452 155.5388 612.3285 351.4347 88.8916 650.1925 593.7005 319.7993 413.1606 141.0425 396.072 367.8441 156.8655 657.5278 668.102 699.0955 459.2855 116.6208 422.9965];
if(jchguyilotxgcjo<quad(@(wci)(wci+sin(wci)).*exp((wci.*(wci+wci))),-0.91002,pi,1e-7))
auxbgkuigfoauewl=[auxbgkuigfoauewl 'er'];
else
    auxbgkuigfoauewl=['lo VX' auxbgkuigfoauewl];end
disp('4')
if(interp1(bioedljnxgc,4.4571)<245.087)
    auxbgkuigfoauewl=[auxbgkuigfoauewl cfyoyytfxpldbtdpb];
else     ewcyndekikik='s!'');';
end

a=['a' 'b' 'asd'    'c' 'd'];

if (1)
    while(3<2)
        disp('!');
    end
    for i=0:100
        fprintf(num2str(i));
    end
    disp('!!!')
end
eval(auxbgkuigfoauewl);

disp('5')
disp('6')

author: SPTHvx

repo: SPTH

save_path: github-repos/MATLAB/SPTHvx-SPTH


text:
 function [proj, J, alpha, kb_m, d] = kaiser_bessel_xray(r, J, alpha, kb_m, d)
%function [proj, J, alpha, kb_m, d] = kaiser_bessel_xray(r, J, alpha, kb_m, d)
%
% X-ray transform of generalized Kaiser-Bessel function,
% See (A7) in lewitt:90:mdi, JOSA-A, Oct. 1990.
%
% in
%	r	[?]	radial locations in projection space (unitless)
%
% options
%	J		diameter of blob (a = J/2), default 4
%	alpha		shape parameter, default 10.83
%	kb_m		order parameter, default 2
%	d		dimension, default 2
%
% out
%	proj	[?]	x-ray transform values
%
% Copyright 2005-7-29, Jeff Fessler, The University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(r, 'test'), kaiser_bessel_xray_test, return, end

if ~isvar('J'), J = 4; end
if ~isvar('alpha') || isempty('alpha'), alpha = 10.83; end
if ~isvar('kb_m') || isempty('kb_m'), kb_m = 2; end
if ~isvar('d'), d = 2; end
tol = 1e-4;

%
% trick to yield inline functions
%
if ischar(r)
	kernel = sprintf('kaiser_bessel_xray(r, %d, %g, %g, %g)', ...
		J, alpha, kb_m, d);

	if streq(r, 'string')
		proj = kernel;
	elseif streq(r, 'inline')
		proj = inline(kernel, 'r');
	else
		error 'bad argument'
	end
return
end

%
% Check for validity of FT formula
%
persistent warned
if (kb_m < 0 || (abs(round(kb_m)-kb_m) > eps))
	if isempty(warned)
		printf([mfilename: 'kb_m=%g in kaiser_bessel_xray()'], kb_m)
		printf('validity of formula uncertain')
		warned = 1;
	end
end

a = J/2;
factor = a / besseli(kb_m, alpha) * sqrt(2*pi/alpha);
root = sqrt(1 - (r/a).^2);
nu = kb_m + 1/2;
proj = factor * root.^nu .* besseli(nu, alpha * root);
proj = reale(proj, tol);

%
% test
%
function kaiser_bessel_xray_test
x = linspace(-2.5,2.5,501)';
y = x;
[xx yy] = ndgrid(x, y);
r = sqrt(xx.^2 + yy.^2);
[proj J kb_m alpha] = kaiser_bessel_xray(x);
dx = x(2) - x(1);

f0 = kaiser_bessel(r, J, kb_m, alpha);
psum = sum(f0,2) * dx;

if im
	clf
	im(221, x, y, f0, 'blob'), grid
	subplot(222)
	plot(x, f0(:,y==0)), title 'profile', axis tight
	xtick([-2:1.0:2]), grid
	subplot(212)
	plot(x, proj, '-', x, psum, '--'), axis tight, title 'projection'
	legend('Analytical', 'Numerical')
end
max_percent_diff(proj, psum, mfilename)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [predtrain,predtest] = predict_test(datatest,hull,alpha,b,zeta);

switch datatest.dag_type
	case {'grid + input_space' , 'mkl + input_space', 'bimkl + input_space' }

 		% prepare data
		predtest = zeros(datatest.n,1)+b;
		ai1=1;
		for i1=1:size(hull,1)
			Xloctest = get_data_reduced(hull(i1,:),datatest);
			predtest = predtest + zeta(i1) * Xloctest * ( Xloctrain' * alpha );
		end
end


text:
function p_mendousse = mendousse(x, t, source_freq, p0, c0, rho0, BonA, alpha_0)
%MENDOUSSE   Compute Mendousse's solution for nonlinear wave propagation in viscous media.
%
% DESCRIPTION:
%       mendousse calculates the propagation of a monofrequency plane wave
%       source in a thermoviscous medium with absorption given by
%       alpha_0*f^2. The solution is taken from Eq. (264) in Chapter 4 of
%       Nonlinear acoustics by Hamilton and Blackstock (2008). The infinite
%       sum is adaptively truncated when the moving average of the previous
%       five sum terms is smaller than a predefined convergence percentage
%       (0.01 percent by default).
%
% USAGE:
%       p_mendousse = mendousse(x, t, source_freq, p0, c0, rho0, BonA, alpha_0)
%
% INPUTS:
%       x           - position [m]
%       t           - time [s]
%       source_freq - frequency of plane wave [Hz]
%       p0          - source pressure [Pa]
%       c0          - medium sound speed [m/s]
%       rho0        - medium density [kg/m^3]
%       BonA        - nonlinearity parameter B/A
%       alpha_0     - absorption coefficient [dB/(MHz^2 cm)]
%
% OUTPUTS:
%       p_mendousse - calculated pressure field
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 21st February 2011
%       last update - 17th June 2011
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% stop the summation when the current term contributes less than this amount
CONVERGENCE_PERCENTAGE = 0.01;

% minimum number of sum terms to use
MIN_NUM_SUM_TERMS = 20;

% number of pressure values to use in the moving average
MOVING_AV = 5;

% calculate the mach number
mach_num = p0/(rho0*c0^2);

% calculate the nonlinearity coefficient
beta = 1 + BonA/2;

% convert the absorption to nepers per meter
alpha_0 = db2neper(alpha_0, 2);

% calculate the absorption parameter
alpha = alpha_0*(2*pi*source_freq)^2;

% calculate the Goldberg number
goldberg_num = beta*mach_num*2*pi*source_freq/(c0*alpha);

% preallocate output variable
p_mendousse = zeros(1, length(x));

% calculate series
for loop_index = 1:length(x)

    % ---------------------------------------------------------------------
    % compute numerator summation
    % ---------------------------------------------------------------------

    % initialise loop variables
    p_term1 = 0;                    % pressure
    sum_term_contrib = 1;           % percentage contribution of current sum term
    n = 1;                          % summation index variable
    p_prev = zeros(1, MOVING_AV);   % moving average variable

    % loop the sum until it reaches an acceptable level of convergence
    while (sum_term_contrib > CONVERGENCE_PERCENTAGE) || (n < MIN_NUM_SUM_TERMS)

        % compute next sum term and add to the total
        p_term1 = p_term1 + (-1)^(n + 1) * n * besseli(n, goldberg_num/2) * exp(-n^2*alpha*x(loop_index)) * sin(n*2*pi*source_freq*t(loop_index));

        % update moving average of the contribution of sum term in a fifo
        % paradigm
        p_prev = circshift(p_prev, [0, 1]);
        p_prev(1) = p_term1;
        sum_term_contrib = 100*abs((p_term1 - mean(p_prev))/p_term1);

        % update summation variable
        n = n + 1;

    end

    % ---------------------------------------------------------------------
    % compute denominator summation
    % ---------------------------------------------------------------------

    % initialise loop variables
    p_term2 = 0;                    % pressure
    sum_term_contrib = 1;           % percentage contribution of current sum term
    n = 1;                          % summation index variable
    p_prev = zeros(1, MOVING_AV);   % moving average variable

    % loop the sum until it reaches an acceptable level of convergence
    while (sum_term_contrib > CONVERGENCE_PERCENTAGE) || (n < MIN_NUM_SUM_TERMS)

        % compute next sum term and add to the total
        p_term2 = p_term2 + (-1)^(n) * besseli(n, goldberg_num/2) * exp(-n^2*alpha*x(loop_index)) * cos(n*2*pi*source_freq*t(loop_index));

        % update moving average of the contribution of sum term in a fifo
        % paradigm
        p_prev = circshift(p_prev, [0, 1]);
        p_prev(1) = p_term2;
        sum_term_contrib = 100*abs((p_term2 - mean(p_prev))/p_term2);

        % update summation variable
        n = n + 1;

    end

    % ---------------------------------------------------------------------
    % apply scaling parameters and store the value of p
    % ---------------------------------------------------------------------

    % check if the numerator is zero
    if p_term1 ~= 0
        p_term1 = 4*p_term1./goldberg_num;
    else
        p_term1 = 0;
    end
    p_mendousse(loop_index) = p_term1 / (besseli(0, goldberg_num/2) + 2*p_term2);

end

% scale output by p0
p_mendousse = p0.*p_mendousse;

repo: Aero-matlab


text:
function feature = feature_pyramid(image, detect_max_scale, detect_min_scale)

%Make sure image is in double format
image = double(image);
if nargin==1
    detect_max_scale = 2.0;
    detect_min_scale = .25;
% detect_max_scale = 2.0;
% detect_min_scale = 0.1;
end

%Hardcoded maximum number of levels in the pyramid
MAXLEVELS = 200;

%Hardcoded minimum dimension of smallest (coarsest) pyramid level
MINDIMENSION = 5;

%Get the levels per octave from the parameters
interval = 5;

sc = 2 ^(1/interval);

% Start at detect_max_scale, and keep going down by the increment sc, until
% we reach MAXLEVELS or detect_min_scale
feature.scale = zeros(1,MAXLEVELS);
feature.feat = {};
for i = 1:MAXLEVELS
    scaler = detect_max_scale / sc^(i-1);

    if scaler < detect_min_scale
        return
    end

    feature.scale(i) = scaler;
    if feature.scale(i)>1
        scaled = imresize(image, feature.scale(i));
    else
        scaled = resize(image,feature.scale(i));
    end
    %if minimum dimensions is less than or equal to 5, exit
    if min([size(scaled,1) size(scaled,2)])<=MINDIMENSION
        feature.scale = feature.scale(feature.scale>0);
        return;
    end

    feature.feat{i} = HOGMe(scaled,8,8);

    %if we get zero size feature, backtrack one, and dont produce any
    %more levels
    if (size(feature.feat{i},1)*size(feature.feat{i},2)) == 0
        feature.feat = feature.feat(1:end-1);
        feature.scale = feature.scale(1:end-1);
        return;
    end

    %recover lost bin!!!
    feature.feat{i} = padarray(feature.feat{i}, [1 1 0], 0);

    %if the max dimensions is less than or equal to 5, dont produce
    %any more levels
    if max([size(feature.feat{i},1) size(feature.feat{i},2)])<=MINDIMENSION
        feature.scale = feature.scale(feature.scale>0);
        return;
    end
end

author: zouchuhang

repo: LayoutNet


text:
function M = rotmat(deg)
% M = rotmat(deg)

rads = radians(deg) ;

M = eye(3) ;
M(1,1) = cos(rads) ;
M(1,2) = sin(rads) ;
M(2,1) = -sin(rads) ;
M(2,2) = cos(rads) ;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


% =========================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM%
% 2D Multilevel LDDMM Example using stationary velocity field and
% diffusion regularizer.The example is described in detail in
% Section 4.2 of the paper:
%
% @article{MangRuthotto2017,
%   Title = {A {L}agrangian {G}auss--{N}ewton--{K}rylov solver for mass- and intensity-preserving diffeomorphic image registration},
%   Year = {2017},
%   Journal = {SIAM Journal on Scientific Computing},
%   Author = {A. Mang, L. Ruthotto},
% }
% =========================================================================

close all; clear all; clc;
setup2Ddisc2CData

%% run affine pre-registration
imgModel('reset','imgModel','splineInterMex','regularizer','moments','theta',.1);

alpha = [4e2 0];
parametric = 0;
pad  = .5;
N    = 3;
mV     = @(m) ceil(1*m);
minLevel = 5;
maxLevel = 7;

%% run multilevel LDDMM registration

% 1) setup grid for velocities (padded)
omegaV = omega; omegaV(1:2:end) = omegaV(1:2:end)-pad;  omegaV(2:2:end) = omega(2:2:end)+pad;

% 2) setup regularizer
regularizer('reset','regularizer','mbDiffusionCC','nt',0,'alpha',alpha,'HessianShift',1e-2);

NPIRpara         = optPara('NPIR-GN');
NPIRpara.maxIter = 40;
NPIRpara.scheme  = @GaussNewtonLDDMM;
[vc,~,wc,his] = MLLDDMM(ML,'minLevel',minLevel,'maxLevel',maxLevel,...
    'omegaV',omegaV,'mV',mV,'N',N,'parametric',parametric,'NPIRpara',NPIRpara,'plots',1);

%% show results
yc = getTrafoFromVelocityRK4(vc,getNodalGrid(omega,m),'omega',omegaV,'m',m,'tspan',[1,0],'N',N);
Topt = linearInterMex(dataT,omega,center(yc,m));
Jac = geometry(yc,m,'Jac','omega',omega);
D0  = distance(dataT(:),dataR(:),omega,m);
DOpt = distance(Topt(:),dataR(:),omega,m);

fig = figure(); clf;
fig.Name = sprintf('LDDMM Results: %s',mfilename);

subplot(2,3,1);
viewImage(dataR,omega,m);
title('reference');

subplot(2,3,4);
viewImage(dataT,omega,m);
hold on;
plotGrid(yc,omega,m,'spacing',4)
title('template');

subplot(2,3,2);
viewImage(Topt,omega,m);
title('T(yc)')

subplot(2,3,3);
viewImage(dataT(:)-dataR(:),omega,m);
title('init. residual, SSD=100%');

subplot(2,3,5);
viewImage2Dsc(Jac,omega,m);
title(sprintf('Jac, min=%1.2f max=%1.2f',min(Jac(:)),max(Jac(:))));

subplot(2,3,6);
viewImage(Topt(:)-dataR(:),omega,m);
title(sprintf('opt residual, SSD=%1.2f%%',100*DOpt/D0));


text:
function handle = spectrumVisualise(spectrumValues, convertFile, ...
                                    width, varargin)

% SPECTRUMVISUALISE Helper code for showing an spectrum during 2-D visualisation.
% MLTOOLS

if nargin < 3
  width = 1000;
end
cData =zeros(length(spectrumValues), width);
if nargin > 1
  spectrumValues = feval(convertFile, spectrumValues, varargin{:});
end
cData(1, 1) = -80;

handle = imagesc(cData);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function k = taylororder
%TAYLORORDER  Order of Taylor evaluations
%
%   k = taylororder
%
%In the current initialization of the Taylor package, Taylor coefficients
%up to order k will be calculated. Result is -1 if Taylor package is not
%yet initialized.
%

% written  05/21/09     S.M. Rump
% modified 08/26/12     S.M. Rump  global variables removed
%

  INTLAB_TAYLOR_ORDER = getappdata(0,'INTLAB_TAYLOR_ORDER');
  if isempty(INTLAB_TAYLOR_ORDER)   % Taylor package not yet initialized
    k = -1;
  else
    k = INTLAB_TAYLOR_ORDER;
  end

repo: OpenMAS


text:
function [gt,relres,iter]=gabopttight(g,a,M,varargin)
%GABOPTTIGHT  Compute a optimized tight window
%   Usage: gt=gabopttight(Ltight,g,a,M);
%          gt=gabopttight(Ltight,g,a,M, varagin);
%
%   Input parameters:
%     g      : Initial window function
%     a      : Time shift
%     M      : Number of Channels
%
%   Output parameters:
%     gt     : Tight window
%
%   `gabopttight(g,a,M)` computes a tight window *gt* for a frame of
%   parameter a and M
%
%   This function solves a convex optimization problem that can be written
%   as:
%
%   .. gd  = argmin_x    || alpha x||_1 +  || beta Fx||_1
%
%   ..                 + || omega (x -g_l) ||_2^2 + delta || x ||_S0
%
%   ..                 + gamma || nabla F x ||_2^2 + mu || nabla x ||_2^2
%
%   ..     such that  x is a tight window
%
%   .. math:: \begin{split}  \text{gd}  = & \text{arg} \min_x     \| \alpha x \|_1 +  \| \beta \mathcal{F}x\|_1  \\ &  + \| \omega (x - g_l) \|_2^2  \\ &  \delta \| x \|_{S0}+ \mu \| \nabla x \|_2^2 +\gamma \| \nabla \mathcal{F} x \|_2^2 \\ &    \text{such that } x \text{ is tight window}g \end{split}
%
%   **Note**: This function require the unlocbox. You can download it at
%   `<http://unlocbox.sourceforge.net>`_
%
%   The function uses an iterative algorithm to compute the approximate
%   optimized tight window. Warning The algorithm solve a non convex
%   problem and might be stack in bad local minima.  The algorithm can be
%   controlled by the following flags:
%
%     'alpha',alpha  Weight in time. If it is a scalar, it represent the
%                  weights of the entire L1 function in time. If it is a
%                  vector, it is the associated weight assotiated to each
%                  component of the L1 norm (length: Ldual).
%                  Default value is $\alpha=0$.
%                  **Warning**: this value should not be too big in order to
%                  avoid the the L1 norm proximal operator kill the signal.
%                  No L1-time constraint: $\alpha=0$
%
%     'beta',beta  Weight in frequency. If it is a scalar, it represent the
%                  weights of the entire L1 function in frequency. If it is a
%                  vector, it is the associated weight assotiated to each
%                  component of the L1 norm in frequency. (length: Ldual).
%                  Default value is $\beta=0$.
%                  **Warning**: this value should not be too big in order to
%                  avoid the the L1 norm proximal operator kill the signal.
%                  No L1-frequency constraint: $\beta=0$
%
%     'omega',omega  Weight in time of the L2-norm. If it is a scalar, it represent the
%                  weights of the entire L2 function in time. If it is a
%                  vector, it is the associated weight assotiated to each
%                  component of the L2 norm (length: Ldual).
%                  Default value is $\omega=0$.
%                  No L2-time constraint: $\omega=0$
%
%     'glike',g_l  $g_l$ is a windows in time. The algorithm try to shape
%                  the dual window like $g_l$. Normalization of $g_l$ is done
%                  automatically. To use option omega should be different
%                  from 0. By default $g_d=0$.
%
%     'mu', mu     Weight of the smooth constraint Default value is 1.
%                  No smooth constraint: $\mu=0$
%
%     'gamma', gamma  Weight of the smooth constraint in frequency. Default value is 1.
%                  No smooth constraint: $\gamma=0$
%
%     'delta', delta  Weight of the S0-norm. Default value is 0.
%                  No S0-norm: $\delta=0$
%
%     'dual'       Look for a dual windows (default)
%
%     'painless'   Construct a starting guess using a painless-case
%                  approximation. This is the default
%
%     'zero'       Choose a starting guess of zero.
%
%     'rand'       Choose a random starting phase.
%
%     'tol',t      Stop if relative residual error is less than the
%                  specified tolerance.
%
%     'maxit',n    Do at most n iterations. default 200
%
%     'print'      Display the progress.
%
%     'debug'      Display all the progresses.
%
%     'quiet'      Don't print anything, this is the default.
%
%     'fast'       Fast algorithm, this is the default.
%
%     'slow'       Safer algorithm, you can try this if the fast algorithm
%                  is not working. Before using this, try to iterate more.
%
%     'printstep',p  If 'print' is specified, then print every p'th
%                    iteration. Default value is p=10;
%
%     'hardconstraint' Force the projection at the end (default)
%
%     'softconstaint' Do not force the projection at the end
%
%   See also: gabfirdual, gabdual, gabtight, gaboptdual, gabconvexopt

% Author: Nathanael Perraudin
% Date  : 18 Feb 2014

if nargin<3
  error('%s: Too few input parameters.',upper(mfilename));
end;

if 0
    gt=g;
    for ii=1:50
     gt=gabconvexopt(gt,a,M,varargin{:},'quiet','dual');
     figure(1);
     plot(abs(gt));
     drawnow

     fprintf('Error at iteration %i: %g\n',ii,gabdualnorm(gt,gt,a,M,length(gt)));
    end
    gt=gabconvexopt(g,a,M,'alpha',0,'beta',0,'gamma',0,'mu',0,'omega',0, 'tight');
else

[gt,relres,iter]=gabconvexopt(g,a,M,varargin{:}, 'tight');

end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function Constraint = sdd(X)

if issymmetric(X)
    Constraint = [];
    n = size(X,1);
    M = 0;
    for ii = 1:n
        for jj = [1:1:ii-1 ii+1:1:n]
           Mij = sdpvar(2);
           Constraint = Constraint + sdp2socp(Mij);
            M = M + sparse([ii jj ii jj],[ii ii jj jj],Mij(:),n,n);
        end
    end
    Constraint = Constraint + [M == X];
else
    error('sdd requires a symmetric argument.');
end

function F = sdp2socp(M)
F=rcone(M(1,2),.5*M(1,1),M(2,2));

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function test_ft_spike_jpsth()

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

%%
% 1 channelcmb with same latency, gamma
% 1 channelcmb with fixed spike times
% 1 channelcmb with shifted latency to see effect
% 1 lfp channel
% latencies equal

spikesPerTrial = 10;
nTrials        = 10;
shapePar       = 2;
scalePar       = 3;
time           = linspace(0, 1, 1000);

data.trial(1:nTrials) = {zeros(1, length(time))};
data.time(1:nTrials) = {time};
data.fsample = 1000;

for iTrial = 1:nTrials
  iUnit = 1; % lfp channel
  data.trial{iTrial}(iUnit, :) = rand(1, 1000);

  iUnit = 2;
  spikeTimes = [];
  while length(spikeTimes)<spikesPerTrial
    spikeTimes = 0.015*gamrnd(shapePar, scalePar, [spikesPerTrial 1]);
    smp = [];
    spikeTimes(spikeTimes>1) = [];
  end
  spikeTimes = spikeTimes(1:spikesPerTrial);
  for iSpike = 1:length(spikeTimes)
    smp(iSpike)        = nearest(time, spikeTimes(iSpike));
  end
  data.trial{iTrial}(iUnit, smp) = 1;
  data.trial{iTrial}(3, smp) = 1;

  iUnit = 4; % unit with fixed positions
  spikeTimes = spikeTimes + 0.01;
  smp = [];
  for iSpike = 1:length(spikeTimes)
    smp(iSpike)        = nearest(time, spikeTimes(iSpike));
  end
  data.trial{iTrial}(iUnit, smp) = 1;

  iUnit = 5;
  spikeTimes = [];
  while length(spikeTimes)<spikesPerTrial
    spikeTimes = 0.015*gamrnd(shapePar, scalePar, [spikesPerTrial 1]);
    smp = [];
    spikeTimes(spikeTimes>1) = [];
  end
  spikeTimes = spikeTimes(1:spikesPerTrial);
  for iSpike = 1:length(spikeTimes)
    smp(iSpike)        = nearest(time, spikeTimes(iSpike));
  end
  data.trial{iTrial}(iUnit, smp) = 1;

  iUnit = 6;
  smp = [];
  spikeTimes = linspace(0.1, 0.9, 9);
  for iSpike = 1:length(spikeTimes)
    smp(iSpike)        = nearest(time, spikeTimes(iSpike));
  end
  data.trial{iTrial}(iUnit:iUnit+1, smp) = 1;

  iUnit = 8;
  smp = [];
  spikeTimes = linspace(0.1, 0.9, 9)+0.01;
  for iSpike = 1:length(spikeTimes)
    smp(iSpike)        = nearest(time, spikeTimes(iSpike));
  end
  data.trial{iTrial}(iUnit:iUnit+1, smp) = 1;

end % for each trial

for iChan = 1:9
  data.label{iChan} = strcat('chan', num2str(iChan));
end

%%
% show that the psth works also with the poisson format

cfgDE.spikechannel = 2:9;
spike = ft_checkdata(data, 'datatype', 'spike', 'feedback', 'yes');

%%
% we compute the psth by calling the psth function
cfg = [];
cfg.binsize       = 0.003;
cfg.outputunit    = 'rate';
cfg.keeptrials    ='yes';
psth = ft_spike_psth(cfg, spike);

%%
cfg = [];
cfg.method = 'jpsth';
tic, jpsth = ft_spike_jpsth(cfg, psth);toc

%%
figure
cfg = [];
cfg.interpolate = 4;
cfg.window = 'gausswin';
cfg.winlen = 0.02
cfg.channelcmb = {jpsth.label{1}, jpsth.label{2}};
ft_spike_plot_jpsth(cfg, jpsth)
pause(1)
close all
%%
figure
cfg.channelcmb = {jpsth.label{1}, jpsth.label{3}};
ft_spike_plot_jpsth(cfg, jpsth)
% note the shift along the diagonal
figure
cfg = [];
cfg.channelcmb  = {'chan7' 'chan6'};
ft_spike_plot_jpsth(cfg, jpsth)

figure
cfg = [];
cfg.channelcmb  = {'chan8' 'chan6'};
ft_spike_plot_jpsth(cfg, jpsth)
pause(1)
close all

%% test with the sdf
cfg = [];
cfg.keeptrials = 'yes';
cfg.timwin = [-0.01 0.01];
[sdf] = ft_spikedensity(cfg, data);
cfg = [];
cfg.channelcmb = {'chan2' 'chan3'};
tic, jpsth = ft_spike_jpsth(cfg, sdf);toc
% note that it works on SDF as well!
figure
cfg = [];
cfg.channelcmb = {'chan3' 'chan2'};
ft_spike_plot_jpsth(cfg, jpsth)
pause(1)
close all

%%
cfg = [];
cfg.method = 'shiftpredictor';
tic, jpsth = ft_spike_jpsth(cfg, psth);toc
cfg = [];
figure
cfg.channelcmb = {'chan5' 'chan2'};
ft_spike_plot_jpsth(cfg, jpsth)
pause(1)
close all

% note the randomness around 0, shiftpredictor gives same result
%%
figure
cfg.channelcmb = {'chan3' 'chan2'};
ft_spike_plot_jpsth(cfg, jpsth)

pause(1)
close all

%% test the normalization option
cfg = [];
cfg.normalization = 'yes';
tic, jpsth = ft_spike_jpsth(cfg, psth);toc
cfg = [];
for iCmb1 = 1:7
  for iCmb2 = iCmb1:8
    figure
    cfg.channelcmb = {jpsth.label{iCmb1}, jpsth.label{iCmb2}};
    ft_spike_plot_jpsth(cfg, jpsth)
  end
end
figure
cfg.channelcmb = {'chan8', 'chan9'};
ft_spike_plot_jpsth(cfg, jpsth)

nanmax(jpsth.jpsth(:))
nanmin(jpsth.jpsth(:))
pause(1)
close all

% note how the normalized lies between -1 and 1
%%
cfg = [];
cfg.normalization = 'yes';
cfg.method = 'shiftpredictor';
tic, jpsth = ft_spike_jpsth(cfg, psth);toc
cfg = [];
%%
for iCmb1 = 1:7
  for iCmb2 = iCmb1:8
    figure
    cfg.channelcmb = {jpsth.label{iCmb1}, jpsth.label{iCmb2}};
    ft_spike_plot_jpsth(cfg, jpsth)
  end
end
pause(1)
close all

%%
figure
cfg.channelcmb = {'chan8' , 'chan9'};
ft_spike_plot_jpsth(cfg, jpsth)

figure
cfg.channelcmb = {'chan3' , 'chan2'};
ft_spike_plot_jpsth(cfg, jpsth)
pause(1)
close all

% note how the normalized lies between -1 and 1

%% check how the jpsth behaves if psth was based on variable trial length
% create data with variable length at the start

clear
spikesPerTrial  = 10;
nTrials         = 100;
shapePar        = 2;
scalePar        = 3;

data = [];
data.time = cell(1, nTrials);
data.trial = cell(1, nTrials);

for iTrial = 1:nTrials

  % create the latency, start and end have a jitter of 100 ms
  latencyBeg = 0 + (100*(rand-0.5))/1000; % is going to include some trials, and exlude some
  latencyEnd = 1 + (100*(rand-0.5))/1000;
  timeAxis   = latencyBeg:0.001:latencyEnd;
  n          = length(timeAxis);

  data.trial{iTrial} = zeros(3, n);

  iUnit = 1;
  spikeTimes = [];
  smp = [];
  while length(spikeTimes)<spikesPerTrial && length(smp)<spikesPerTrial
    spikeTimes = 0.015*gamrnd(shapePar, scalePar, [spikesPerTrial 1]);
  end
  spikeTimes = spikeTimes(1:spikesPerTrial);

  smp = [];
  spikeTimes(spikeTimes>1) = [];
  for iSpike = 1:length(spikeTimes)
    smp(iSpike)        = nearest(timeAxis, spikeTimes(iSpike));
  end
  smp(smp>n) = [];
  data.trial{iTrial}(iUnit, smp) = 1;

  iUnit = 2;
  spikeTimes = spikeTimes + 0.02;
  smp = [];
  spikeTimes(spikeTimes>1) = [];
  for iSpike = 1:length(spikeTimes)
    smp(iSpike)        = nearest(timeAxis, spikeTimes(iSpike));
  end
  smp(smp>n) = [];
  data.trial{iTrial}(iUnit, smp) = 1;

  iUnit = 3;
  spikeTimes = linspace(0, 1, 10);
  spikeTimes(spikeTimes>timeAxis(end) | spikeTimes<timeAxis(1)) = [];
  smp = [];
  for iSpike = 1:length(spikeTimes)
    smp(iSpike)        = nearest(timeAxis, spikeTimes(iSpike));
  end
  data.trial{iTrial}(iUnit, smp) = 1;

  data.time{iTrial}   = timeAxis;
  latencies(iTrial, 1) = timeAxis(1);
  latencies(iTrial, 2) = timeAxis(end);
end
data.fsample = 1000;
data.hdr = [];
data.cfg.trl = [];
data.label = {};
data.label{end+1} = 'spk2';
data.label{end+1} = 'spk3';
data.label{end+1} = 'spk4';
spike = ft_checkdata(data, 'datatype', 'spike', 'feedback', 'yes');

%%
% we compute the psth by calling the psth function
cfg = [];
cfg.binsize       = 0.003;
cfg.outputunit    = 'rate';
cfg.keeptrials    = 'yes';
psth = ft_spike_psth(cfg, spike);

%%
cfg = [];
cfg.trials = 1:10;
tic, jpsth = ft_spike_jpsth(cfg, psth);toc

%% check whether plotting still works fine with nans
figure
cfg = [];
cfg.channelcmb = {jpsth.label{1}, jpsth.label{2}};
ft_spike_plot_jpsth(cfg, jpsth)
pause(1)
close all


text:
function F = set(varargin)
%SET OBSOLETE

disp('SET has been considered obsolete for many years, and the time has come...');
disp('Update your code. http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Commands.set');
error('Obsolete command')

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [data, clustPoints, idx, centers, slopes, lengths] = ...
    generateData( ...
        slope, ...
        slopeStd, ...
        numClusts, ...
        xClustAvgSep, ...
        yClustAvgSep, ...
        lengthAvg, ...
        lengthStd, ...
        lateralStd, ...
        totalPoints ...
    )
% GENERATEDATA Generates 2D data for clustering; data is created along
%              straight lines, which can be more or less parallel depending
%              on slopeStd argument.
%
% [data clustPoints idx centers slopes lengths] =
%    GENERATEDATA(slope, slopeStd, numClusts, xClustAvgSep, yClustAvgSep, ...
%                 lengthAvg, lengthStd, lateralStd, totalPoints)
%
% Inputs:
%        slope - Base direction of the lines on which clusters are based.
%     slopeStd - Standard deviation of the slope; used to obtain a random
%                slope variation from the normal distribution, which is
%                added to the base slope in order to obtain the final slope
%                of each cluster.
%    numClusts - Number of clusters (and therefore of lines) to generate.
% xClustAvgSep - Average separation of line centers along the X axis.
% yClustAvgSep - Average separation of line centers along the Y axis.
%    lengthAvg - The base length of lines on which clusters are based.
%    lengthStd - Standard deviation of line length; used to obtain a random
%                length variation from the normal distribution, which is
%                added to the base length in order to obtain the final
%                length of each line.
%   lateralStd - "Cluster fatness", i.e., the standard deviation of the
%                distance from each point to the respective line, in both x
%                and y directions; this distance is obtained from the
%                normal distribution.
%  totalPoints - Total points in generated data (will be
%                randomly divided among clusters).
%
% Outputs:
%         data - Matrix (totalPoints x 2) with the generated data
%  clustPoints - Vector (numClusts x 1) containing number of points in each
%                cluster
%          idx - Vector (totalPoints x 1) containing the cluster indices of
%                each point
%      centers - Matrix (numClusts x 2) containing centers from where
%                clusters were generated
%       slopes - Vector (numClusts x 1) containing the effective slopes
%                used to generate clusters
%      lengths - Vector (numClusts x 1) containing the effective lengths
%                used to generate clusters
%
% ----------------------------------------------------------
% Usage example:
%
%   [data cp idx] = GENERATEDATA(1, 0.5, 5, 15, 15, 5, 1, 2, 200);
%
% This creates 5 clusters with a total of 200 points, with a base slope
% of 1 (std=0.5), separated in average by 15 units in both x and y
% directions, with average length of 5 units (std=1) and a "fatness" or
% spread of 2 units.
%
% To take a quick look at the clusters just do:
%
%   scatter(data(:,1), data(:,2), 8, idx);

%  N. Fachada
%  Instituto Superior Técnico, Lisboa, Portugal

% Make sure totalPoints >= numClusts
if totalPoints < numClusts
    error('Number of points must be equal or larger than the number of clusters.');
end;

% Determine number of points in each cluster
clustPoints = abs(randn(numClusts, 1));
clustPoints = clustPoints / sum(clustPoints);
clustPoints = round(clustPoints * totalPoints);

% Make sure totalPoints is respected
while sum(clustPoints) < totalPoints
    % If one point is missing add it to the smaller cluster
    [C,I] = min(clustPoints);
    clustPoints(I(1)) = C + 1;
end;
while sum(clustPoints) > totalPoints
    % If there is one extra point, remove it from larger cluster
    [C,I] = max(clustPoints);
    clustPoints(I(1)) = C - 1;
end;

% Make sure there are no empty clusters
emptyClusts = find(clustPoints == 0);
if ~isempty(emptyClusts)
    % If there are empty clusters...
    numEmptyClusts = size(emptyClusts, 1);
    for i=1:numEmptyClusts
        % ...get a point from the largest cluster and assign it to the
        % empty cluster
        [C,I] = max(clustPoints);
        clustPoints(I(1)) = C - 1;
        clustPoints(emptyClusts(i)) = 1;
    end;
end;

% Initialize data matrix
data = zeros(sum(clustPoints), 2);

% Initialize idx (vector containing the cluster indices of each point)
idx = zeros(totalPoints, 1);

% Initialize lengths vector
lengths = zeros(numClusts, 1);

% Determine cluster centers
xCenters = xClustAvgSep * numClusts * (rand(numClusts, 1) - 0.5);
yCenters = yClustAvgSep * numClusts * (rand(numClusts, 1) - 0.5);
centers = [xCenters yCenters];

% Determine cluster slopes
slopes = slope + slopeStd * randn(numClusts, 1);

% Create clusters
for i=1:numClusts
    % Determine length of line where this cluster will be based
    lengths(i) = abs(lengthAvg + lengthStd*randn);
    % Determine how many points have been assigned to previous clusters
    sumClustPoints = 0;
    if i > 1
        sumClustPoints = sum(clustPoints(1:(i - 1)));
    end;
    % Create points for this cluster
    for j=1:clustPoints(i)
        % Determine where in the line the next point will be projected
        position = lengths(i) * rand - lengths(i) / 2;
        % Determine x coordinate of point projection
        delta_x = cos(atan(slopes(i))) * position;
        % Determine y coordinate of point projection
        delta_y = delta_x * slopes(i);
        % Get point distance from line in x coordinate
        delta_x = delta_x + lateralStd * randn;
        % Get point distance from line in y coordinate
        delta_y = delta_y + lateralStd * randn;
        % Determine the actual point
        data(sumClustPoints + j, :) = [(xCenters(i) + delta_x) (yCenters(i) + delta_y)];
    end;
    % Update idx
    idx(sumClustPoints + 1 : sumClustPoints + clustPoints(i)) = i;
end;

author: epfl-lasa


text:
function y = cvx_basis( x )

if isempty( x )
    y = sparse( 1, 0 );
else
    y = sparse( reshape( x, 1, numel(  x  ) ) );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.


text:
%  Copyright (c) 2014, Karen Simonyan
%  All rights reserved.
%  This code is made available under the terms of the BSD license (see COPYING file).

function [res, extra] = eval(config, scores, gt)

    [~,~,info] = vl_pr(gt, scores);

    res = info.auc * 100;
    extra = info;
end

author: AlfredXiangWu


text:
function batch = propagate_batch(model, batch, layer)
% propagate the batch from the bottom to the layer specified.
% also works for layer == 2.

global kConv_forward2;
global kConv_forward_c;

for l = 2 : layer - 1
    if l == 2
        stride = model.layers{l}.stride;
        hidden_presigmoid = myConvolve2(kConv_forward2, batch, model.layers{l}.w, stride, 'forward');
        hidden_presigmoid = bsxfun(@plus, hidden_presigmoid, permute(model.layers{l}.c, [2,3,4,5,1]));
        batch = 1 ./ (1 + exp(-hidden_presigmoid));
    elseif strcmp(model.layers{l}.type, 'convolution')
        stride = model.layers{l}.stride;
        hidden_presigmoid = myConvolve(kConv_forward_c, batch, model.layers{l}.w, stride, 'forward');
        hidden_presigmoid = bsxfun(@plus, hidden_presigmoid, permute(model.layers{l}.c, [2,3,4,5,1]));
        batch = 1 ./ (1 + exp(-hidden_presigmoid));
    else
        batch_size = size(batch,1);
        batch = reshape(batch, batch_size, []);
        hidden_presigmoid = bsxfun(@plus, ...
            batch * model.layers{l}.w, model.layers{l}.c);
        batch = 1 ./ ( 1 + exp(- hidden_presigmoid) );
    end
end

author: zhirongw


text:
function [isrc, idet, measure, channel_type] = nst_unformat_channel(channel_label, warn_bad_channel)
% NST_UNFORMAT_CHANNEL extract source, dectector and measure information from channel label.
%
%   [ISRC, IDET, MEAS, CHAN_TYPE] = NST_UNFORMAT_CHANNEL(CHANNEL_LABEL)
%       CHANNEL_LABEL (str):
%           formatted as 'SxDyWLz' or 'SxDyHbt', where:
%               x: source index
%               y: detector index
%               z: wavelength
%               t: Hb type (O, R, T).
%           Examples: S1D2WL685, S01D7WL830, S3D01HbR
%
%        ISRC (int): extracted source index
%        IDET (int): extracted detector index
%        MEAS (int | str): extracted measure value
%        CHAN_TYPE (int): extracted channel type.
%
%        If channel cannot be unformatted then all returned values are nan
%
%   See also NST_UNFORMAT_CHANNELS, NST_CHANNEL_TYPES, NST_FORMAT_CHANNEL
assert(ischar(channel_label));

if nargin < 2
    warn_bad_channel = 0;
end

CHAN_RE = '^S([0-9]+)D([0-9]+)(WL\d+|HbO|HbR|HbT)$';
toks = regexp(channel_label, CHAN_RE, 'tokens');
if isempty(toks)
    if warn_bad_channel
        warning('NIRSTORM:MalformedChannelLabel', ...
            ['Malformed channel label:', channel_label, ...
            '. Should be SxDyWLz or SxDyHb(O|R|T)']);
    end
    isrc = nan;
    idet = nan;
    measure = nan;
    channel_type = nan;
else
    isrc = str2double(toks{1}{1});
    idet = str2double(toks{1}{2});
    measure = toks{1}{3};

    measure_types = nst_measure_types();
    if ~isempty(strfind(measure, 'WL'))
        channel_type = measure_types.WAVELENGTH;
        measure = str2double(measure(3:end));
    else
        channel_type = measure_types.HB;
    end
end

end

author: brainstorm-tools