text:
% STD_DIPPLOT - Commandline function to plot cluster component dipoles. Dipoles for each
%                 named cluster is displayed in a separate figure. To view all the clustered
%                 components in the STUDY on the same figure (in a separate subplot), all
%                 STUDY clusters must be requested.
%                 To visualize dipoles, they first must be stored in the EEG dataset structures
%                 using DIPFIT. Only components that have dipole locations will be displayed,
%                 along with the cluster mean dipole (in red).
% Usage:
%                 >> [STUDY] = std_dipplot(STUDY, ALLEEG, clusters);
% Inputs:
%   STUDY      - EEGLAB STUDY set comprising some or all of the EEG datasets in ALLEEG.
%   ALLEEG     - global EEGLAB vector of EEG structures for the dataset(s) included in
%                the STUDY. ALLEEG for a STUDY set is typically created using LOAD_ALLEEG.
%
% Optional inputs:
%   'clusters' - [numeric vector | 'all']  -> specific cluster numbers to plot.
%                'all'  -> plot all clusters in STUDY.
%                {default: 'all'}.
%   'comps'    - [numeric vector]  -> indices of the cluster components to plot.
%                'all'  -> plot all the components in the cluster
%                {default: 'all'}.
%   'mode'     - ['together'|'apart'|'multicolor'] Display all requested cluster on one
%                figure ('together') or separate figures ('apart').
%                'together'-> plot all 'clusters' individually in one multi-pane figure (without the gui).
%                'apart'   -> plot each cluster in a separate figure.
%                'multicolor' -> plot all clusters in one figure,
%                Note that this parameter has no effect if the 'comps' option (above) is used.
%                {default: 'together'}
%   'figure'   - ['on'|'off'] plots on a new figure ('on')  or plots on current
%                figure ('off'). If 'figure','off' does not display gui controls,
%                Useful for incomporating a cluster dipplot into a complex figure.
%                {default: 'on'}.
%   'groups'   - ['on'|'off'] use different colors for different groups.
%                {default: 'off'}.
%   'dipcolor' - [cell vector] color for dipoles in each cluster. (multicolor mode)
%   'dipsize'  - [numeric vector] size for each cluster. (multicolor mode)
%                {default if unspecified: will automatically color/size each
%                cluster}
% Outputs:
%   STUDY      - the input STUDY set structure modified with plotted cluster
%                mean dipole, to allow quick replotting (unless cluster means
%                already exists in the STUDY).
%   Example:
%   >> [STUDY] = std_dipplot(STUDY,ALLEEG, 'clusters', 5, 'mode', 'apart', 'figure', 'off');
%                % Plot cluster-5 component dipoles (in blue), plus their mean dipole (in red),
%                % on an existing (gui-less) figure.
%
%  See also  POP_CLUSTEDIT, DIPPLOT
%
% Authors:  Hilit Serby, Arnaud Delorme, Scott Makeig, SCCN, INC, UCSD, June, 2005
%          'groups' added by Makoto Miyakoshi on June 2012.
%          'multicolor' mode added by John Iversen to draw all clusters on
%           a single panel, with each cluster indicated by different color/size dipoles.

% Copyright (C) Hilit Serby, SCCN, INC, UCSD, June 08, 2005, hilit@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function STUDY = std_dipplot(STUDY, ALLEEG, varargin)

% Set default values
cls = []; % plot all clusters in STUDY
figureon = 1; % plot on a new figure
mode = 'apart';

STUDY = pop_dipparams(STUDY, 'default');

opt_dipplot = {'projlines',STUDY.etc.dipparams.projlines, 'axistight', STUDY.etc.dipparams.axistight, 'projimg', STUDY.etc.dipparams.projimg, 'dipolelength', 0, 'density', STUDY.etc.dipparams.density};

dipcolor = [];
dipsize = [];
%, 'spheres', 'on'
groupval = 'off';
nosphere = true;
for k = 3:2:nargin
    switch varargin{k-2}
        case 'clusters'
            if isnumeric(varargin{k-1})
                cls = varargin{k-1};
                if isempty(cls)
                    cls = 2:length(STUDY.cluster);
                end
            else
                if ischar(varargin{k-1}) && strcmpi(varargin{k-1}, 'all')
                    cls = 2:length(STUDY.cluster);
                else
                    error('std_dipplot: ''clusters'' input takes either specific clusters (numeric vector) or keyword ''all''.');
                end
            end
            if length(cls) == 1, mode = 'apart'; else mode = 'together'; end
        case 'comps'
            if strcmpi(STUDY.etc.dipparams.density, 'on')
                disp('Single dipole should not be plotted using dipole density, reverting to dipole plotting');
                opt_dipplot{end} = 'off';
            end
            STUDY = std_plotcompdip(STUDY, ALLEEG,  cls, varargin{k-1}, opt_dipplot{:});
            return;
        case 'plotsubjects', % do nothing
        case 'mode', mode = varargin{k-1};
        case 'groups', groupval = varargin{k-1};
        case 'figure'
            if strcmpi(varargin{k-1},'off')
                if strcmpi(STUDY.etc.dipparams.density, 'on')
                    disp('Cannot plot dipole density within figure, reverting to dipole plotting');
                    opt_dipplot{end} = 'off';
                end
                opt_dipplot{end + 1} = 'gui';
                opt_dipplot{end + 1} = 'off';
                figureon = 0;
            end
      case 'dipcolor'
        dipcolor = varargin{k-1};
      case 'dipsize'
        dipsize = varargin{k-1};
      case 'spheres'
        opt_dipplot = { opt_dipplot{:}, 'spheres',  varargin{k-1} };
        nosphere = false;
    end
end
if nosphere
    opt_dipplot = { opt_dipplot{:}, 'spheres', 'on' };
end

if strcmpi(mode, 'together')
    if strcmpi(STUDY.etc.dipparams.density, 'on')
        disp('Cannot plot dipole density within figure, reverting to dipole plotting');
        opt_dipplot{end} = 'off';
    end
end

% select clusters to plot
% -----------------------
if isempty(cls)
    tmp =[];
    cls = 2:length(STUDY.cluster); % plot all clusters in STUDY
    for k = 1: length(cls)
        % don't include 'Notclust' clusters
        if ~strncmpi('Notclust',STUDY.cluster(cls(k)).name,8) && ~strncmpi('ParentCluster',STUDY.cluster(cls(k)).name,13)
            tmp = [tmp cls(k)];
        end
    end
    cls = tmp;
end

if strcmpi(mode, 'apart')  % case each cluster on a separate figure
    for clus = 1: length(cls) % For each cluster requested
        if length(STUDY.cluster(cls(clus)).comps) > 0  % check there are comps in cluster
            max_r = 0;
            clear cluster_dip_models;
            len = length(STUDY.cluster(cls(clus)).comps);
            ndip = 0;
            dip_ind = [];
            if ~isfield(STUDY.cluster(cls(clus)),'dipole')
                STUDY = std_centroid(STUDY,ALLEEG, cls(clus) , 'dipole');
            elseif isempty(STUDY.cluster(cls(clus)).dipole)
                STUDY = std_centroid(STUDY,ALLEEG, cls(clus) , 'dipole');
            end
            for k = 1:len
                abset   = STUDY.datasetinfo(STUDY.cluster(cls(clus)).sets(1,k)).index;
                subject = STUDY.datasetinfo(STUDY.cluster(cls(clus)).sets(1,k)).subject;
               if ~isfield(ALLEEG(abset), 'dipfit')
                   warndlg2(['No dipole information available in dataset ' ALLEEG(abset).filename ' , abort plotting'], 'Aborting plot dipoles');
                   return;
               end
               comp = STUDY.cluster(cls(clus)).comps(k);
               cluster_dip_models(k).posxyz = ALLEEG(abset).dipfit.model(comp).posxyz;
               cluster_dip_models(k).momxyz = ALLEEG(abset).dipfit.model(comp).momxyz;
               cluster_dip_models(k).rv = ALLEEG(abset).dipfit.model(comp).rv;
               if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
                   if isfield(ALLEEG(abset).dipfit, 'hdmfile') %dipfit 2 spherical model
                       load('-mat', ALLEEG(abset).dipfit.hdmfile);
                       max_r = max(max_r, max(vol.r));
                   else % old version of dipfit
                       max_r = max(max_r,max(ALLEEG(abset).dipfit.vol.r));
                   end
               end
               comp_to_disp{k} = [subject  ', ' 'IC' num2str(comp) ];
               if ~isempty(cluster_dip_models(k).posxyz)
                   ndip = ndip +1;
                   dip_ind = [dip_ind k];
               end
            end % finished going over cluster comps

            STUDY.cluster(cls(clus)).dipole = computecentroid(cluster_dip_models);
            cluster_dip_models(end + 1) = STUDY.cluster(cls(clus)).dipole;

           % additional options
           % ------------------
           dip_color = cell(1,ndip+1);
           dip_color(1:ndip) = {'b'};
           dip_color(end) = {'r'};
           options = opt_dipplot;
           options{end+1} =  'mri';
           options{end+1} =  ALLEEG(abset).dipfit.mrifile;
           options{end+1} =  'coordformat';
           options{end+1} =  ALLEEG(abset).dipfit.coordformat;
           options{end+1} =  'dipnames';
           options{end+1} = {comp_to_disp{dip_ind } [STUDY.cluster(cls(clus)).name ' mean']};
           options{end+1} = 'color';
           options{end+1} = dip_color;

           % if 'groups'==1, overwrite cluster_dip_models, dip_color and dipnames in option -makoto
           if strcmpi(groupval, 'on')
               [cluster_dip_models, options] = dipgroups(ALLEEG, STUDY, cls, comp_to_disp, cluster_dip_models, options);
               break
           end

           if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
               options{end+1} = 'sphere';
               options{end+1} = max_r;
           else
               options{end+1} = 'meshdata';
               options{end+1} = ALLEEG(abset).dipfit.hdmfile;
           end
           %enable both lines and images to be projected; increase #dipole limits
           if ndip < 20 && strcmpi(options{1}, 'projlines') && length(cls) == 1
               options{2} = 'on';
           end
           if ndip < 20 && strcmpi(options{5}, 'projimg') && length(cls) == 1
             options{6} = 'on';
           end

           % Dealing with projection lines
           if strcmpi(options{2},'on')
               projlinesvect = ones(1,length(cluster_dip_models));
           elseif strcmpi(options{2},'off') && strcmp(STUDY.etc.dipparams.centrline,'on')
               projlinesvect = zeros(1,length(cluster_dip_models));
               projlinesvect(end) = 1;
           elseif strcmpi(options{2},'off')
               projlinesvect = zeros(1,length(cluster_dip_models));
           end
           options{2} = projlinesvect;

           if figureon
               dipplot(cluster_dip_models, options{:});
               fig_h = gcf;
               set(fig_h,'Name', [STUDY.cluster(cls(clus)).name ' - ' num2str(length(unique(STUDY.cluster(cls(clus)).sets(1,:)))) ...
                       ' sets - ' num2str(length(STUDY.cluster(cls(clus)).comps)) ' components (' num2str(ndip) ' dipoles)' ],'NumberTitle','off');
           else
               dipplot(cluster_dip_models, options{:},'view', [0.5 -0.5 0.5]);
               for gind = 1:length(options) % remove the 'gui' 'off' option
                   if ischar(options{gind})
                       if strfind(options{gind}, 'gui')
                           break;
                       end
                   end
               end
               options(gind:gind+1) = [];
               dipinfo.dipmod =  cluster_dip_models;
               dipinfo.op = options;
               diptitle = [STUDY.cluster(cls(clus)).name ', ' num2str(length(unique(STUDY.cluster(cls(clus)).sets(1,:)))) ' sets -' ...
                   num2str(length(STUDY.cluster(cls(clus)).comps)) ' components (' num2str(ndip) ' dipoles)' ];
               dipinfo.title = diptitle;
               set(gcf, 'UserData', dipinfo);
               set(gca,'UserData', dipinfo);
               rotate3d off;
               axcopy(gca, ['dipinfo = get(gca, ''''UserData''''); dipplot(dipinfo.dipmod, dipinfo.op{:}); set(gcf, ''''Name'''', dipinfo.title,''''NumberTitle'''',''''off''''); ']);
           end
        end % finished the if condition that cluster isn't empty
    end % finished going over requested clusters
end

if strcmpi(mode, 'together')  % case all clusters are plotted in the same figure (must be a new figure)
    N = length(cls);
    rowcols(2) = ceil(sqrt(N)); % Number of rows in the subplot figure.
    rowcols(1) = ceil(N/rowcols(2));
    fig_h = figure;
    orient tall
    set(fig_h,'Color', 'black');
    set(fig_h,'Name', 'All clusters dipoles','NumberTitle','off');
    set(fig_h, 'resize','off');
    for l = 1:N
        len = length(STUDY.cluster(cls(l)).comps);
        max_r = 0;
        clear cluster_dip_models;
        if ~isfield(STUDY.cluster(cls(l)),'dipole')
            STUDY = std_centroid(STUDY,ALLEEG, cls(l), 'dipole');
        elseif isempty(STUDY.cluster(cls(l)).dipole)
            STUDY = std_centroid(STUDY,ALLEEG, cls(l), 'dipole');
        end
        for k = 1: len
            abset = STUDY.datasetinfo(STUDY.cluster(cls(l)).sets(1,k)).index;
           if ~isfield(ALLEEG(abset), 'dipfit')
               warndlg2(['No dipole information available in dataset ' num2str(abset) ' , abort plotting'], 'Aborting plot dipoles');
               return;
           end
           comp = STUDY.cluster(cls(l)).comps(k);
           cluster_dip_models(k).posxyz = ALLEEG(abset).dipfit.model(comp).posxyz;
           cluster_dip_models(k).momxyz = ALLEEG(abset).dipfit.model(comp).momxyz;
           cluster_dip_models(k).rv = ALLEEG(abset).dipfit.model(comp).rv;
           if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
                if isfield(ALLEEG(abset).dipfit, 'hdmfile') %dipfit 2 spherical model
                   load('-mat', ALLEEG(abset).dipfit.hdmfile);
                   max_r = max(max_r, max(vol.r));
               else % old version of dipfit
                   max_r = max(max_r,max(ALLEEG(abset).dipfit.vol.r));
               end
           end
        end % finished going over cluster comps
        STUDY.cluster(cls(l)).dipole = computecentroid(cluster_dip_models);
        cluster_dip_models(end + 1) = STUDY.cluster(cls(l)).dipole;
        dip_color = cell(1,length(cluster_dip_models));
        dip_color(1:end-1) = {'b'};
        dip_color(end) = {'r'};
        options = opt_dipplot;
        options{end + 1} =  'gui';
        options{end + 1} =  'off';
        options{end+1} =  'mri';
        options{end+1} =  ALLEEG(abset).dipfit.mrifile;
        options{end+1} =  'coordformat';
        options{end+1} =  ALLEEG(abset).dipfit.coordformat;
        options{end+1} = 'color';
        options{end+1} = dip_color;
        if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
            options{end+1} = 'sphere';
            options{end+1} = max_r;
        else
            options{end+1} = 'meshdata';
            options{end+1} = ALLEEG(abset).dipfit.hdmfile;
        end

        % Dealing with projection lines
        if strcmpi(options{2},'on')
            projlinesvect = ones(1,length(cluster_dip_models));
        elseif strcmpi(options{2},'off') && strcmp(STUDY.etc.dipparams.centrline,'on')
            projlinesvect = zeros(1,length(cluster_dip_models));
            projlinesvect(end) = 1;
        elseif strcmpi(options{2},'off')
            projlinesvect = zeros(1,length(cluster_dip_models));
        end
           options{2} = projlinesvect;

        hsbplot = subplot(rowcols(1),rowcols(2),l);
        % Creating new axis for title only
        set(hsbplot,'visible','off');
        htitle = title([ STUDY.cluster(cls(l)).name ' (' num2str(length(unique(STUDY.cluster(cls(l)).sets(1,:)))) ' Ss, ' num2str(length(STUDY.cluster(cls(l)).comps)),' ICs)'],'color','white','Visible', 'on');
        if rowcols(1)> 1, set(htitle,'Position', get(htitle,'Position')+ [0 0.05 0].*get(htitle,'Position') ); end
        axes('Position', get(hsbplot,'Position'),'Color', 'black');
        dipplot(cluster_dip_models, options{:});

        %diptitle = [STUDY.cluster(cls(l)).name ', ' num2str(length(unique(STUDY.cluster(cls(l)).sets(1,:)))) 'Ss'];
        %title(diptitle, 'Color', 'white');
        % Complex axcopy
        %if l == 1
        %    for gind = 1:length(options) % remove the 'gui' 'off' option
        %        if ischar(options{gind})
        %            if strfind(options{gind}, 'gui')
        %                break;
        %            end
        %        end
        %    end
        %    options(gind:gind+1) = [];
        %end
        %dipinfo.dipmod =  cluster_dip_models;
        %dipinfo.op = options;
        %dipinfo.title = diptitle;
        %set(gcf, 'UserData', dipinfo);
        %set(gca,'UserData', dipinfo);
        %axcopy(gcf, ['dipinfo = get(gca, ''''UserData''''); dipplot(dipinfo.dipmod, dipinfo.op{:}); set(gcf, ''''Name'''', dipinfo.title,''''NumberTitle'''',''''off'''');']);
   end %finished going over all clusters
   set(fig_h, 'resize','on');
end % finished case of 'all' clusters

% ========================================================================================
% multicolor mode
%   all clusters are plotted in the same axis, with each cluster indicated by color/size
%   also enable 'data cursor' so that in data tip mode, clicking on a dipole
%    will display its name
if strcmpi(mode, 'multicolor')
  N = length(cls);
  %%%%%%%%%%%%%%%%%%%%% color list %%%%%%%%%%%%%%%%%%%%%
  % This color list was developed for std_envtopo
  % modified from dipgroups below
  colors{1}  = [1 1 1];            % White
  colors{2}  = [1 1 0];            % Yellow
  colors{3}  = [1 0 1];            % Fuchsia
  colors{4}  = [1 0 0];            % Red
  colors{5}  = [0.875 0.875 0.875]; % Silver
  colors{6}  = [0.5 0.5 0.5];      % Gray
  colors{7}  = [0.5 0.5 0];        % Olive
  colors{8}  = [0.5 0 0.5];        % Purple
  colors{9}  = [0.5 0 0];          % Maroon
  colors{10} = [0 1 1];            % Aqua
  colors{11} = [0 1 0];            % Lime
  colors{12} = [0 0.5 0.5];        % Teal
  colors{13} = [0 0.5 0];          % Green
  colors{14} = [0 0 1];            % Blue
  colors{15} = [0 0 0.5];          % Navy
  colors{16} = [0 0 0];            % Black
  % Choosing and sorting 13 colors for clusters: Red, Green, Blue,
  % Fuchsia, Lime, Aqua, Maroon, Olive, Purple, Teal, Navy, Gray, and White
  colors = colors([4 13 14 3 11 10 9 7 8 12 15 6 1 ]);
  fig_h = figure;
  orient tall
  set(fig_h,'Color', 'black');
  set(fig_h,'Name', 'All clusters dipoles','NumberTitle','off');
  set(fig_h, 'resize','off');

  idx = 0; %cumulative dipole index
  centroidIdx = [];
  clear cluster_dip_models;
  for l = 1:N %loop over clusters
    len = length(STUDY.cluster(cls(l)).comps);
    max_r = 0;

    %color, size for every cluster can be passed as arguments, or
    %automatically iterate through a set of colors/sizes
    if ~isempty(dipcolor)
      clusterColors{l} = dipcolor{idx};
    else
      colorIndex = mod(l-1, length(colors))+1;
      clusterColors{l} = colors{colorIndex};
    end

    if ~isempty(dipsize)
      clusterSizes(l) = dipsize(idx);
    else
      %after rotating through color list once, change size of dipole
      if l <= length(colors)
        clusterSizes(l) = 30;
      elseif l <= 2*length(colors)
        clusterSizes(l) = 20;
      else
        clusterSizes(l) = 15;
      end
    end

    if ~isfield(STUDY.cluster(cls(l)),'dipole')
      STUDY = std_centroid(STUDY,ALLEEG, cls(l), 'dipole');
    elseif isempty(STUDY.cluster(cls(l)).dipole)
      STUDY = std_centroid(STUDY,ALLEEG, cls(l), 'dipole');
    end
    clustStartIdx = idx + 1;

    for k = 1:len %loop over components within a cluster
      idx = idx + 1;
      abset = STUDY.datasetinfo(STUDY.cluster(cls(l)).sets(1,k)).index;
      subjname = STUDY.datasetinfo(STUDY.cluster(cls(l)).sets(1,k)).subject;
      if ~isfield(ALLEEG(abset), 'dipfit')
        warndlg2(['No dipole information available in dataset ' num2str(abset) ' , abort plotting'], 'Aborting plot dipoles');
        return;
      end
      comp = STUDY.cluster(cls(l)).comps(k);
      cluster_dip_models(idx).posxyz = ALLEEG(abset).dipfit.model(comp).posxyz;
      cluster_dip_models(idx).momxyz = ALLEEG(abset).dipfit.model(comp).momxyz;
      cluster_dip_models(idx).rv = ALLEEG(abset).dipfit.model(comp).rv;
      if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
        if isfield(ALLEEG(abset).dipfit, 'hdmfile') %dipfit 2 spherical model
          load('-mat', ALLEEG(abset).dipfit.hdmfile);
          max_r = max(max_r, max(vol.r));
        else % old version of dipfit
          max_r = max(max_r,max(ALLEEG(abset).dipfit.vol.r));
        end
      end

      dip_color{idx} = clusterColors{l};
      dip_size(idx) = clusterSizes(l);
      dip_label{idx} = sprintf('Cls %d (%s IC%d)',cls(l), subjname, comp);
    end % finished going over cluster comps

    %add the cluster centroid
    clustEndIdx = idx;
    STUDY.cluster(cls(l)).dipole = computecentroid(cluster_dip_models(clustStartIdx:clustEndIdx));
    idx = idx + 1;
    centroidIdx(end+1) = idx;
    cluster_dip_models(idx) = STUDY.cluster(cls(l)).dipole;
    dip_color(idx) = {'k'};
    dip_size(idx) = 10;
    dip_label{idx} = sprintf('Cls %d centroid',cls(l));
    clusterLabels{l} = sprintf('Cls %d',cls(l));

  end %loop over clusters

  options = opt_dipplot;
  options{end + 1} =  'gui';
  options{end + 1} =  'off';
  options{end+1} =  'mri';
  options{end+1} =  ALLEEG(abset).dipfit.mrifile;
  options{end+1} =  'coordformat';
  options{end+1} =  ALLEEG(abset).dipfit.coordformat;
  options{end+1} = 'color';
  options{end+1} = dip_color;
  options{end+1} = 'dipolesize';
  options{end+1} = dip_size;
  options{end+1} = 'dipnames';
  options{end+1} = dip_label;
  if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
    options{end+1} = 'sphere';
    options{end+1} = max_r;
  else
    options{end+1} = 'meshdata';
    options{end+1} = ALLEEG(abset).dipfit.hdmfile;
  end

  % Dealing with projection lines
  if strcmpi(options{2},'on')
      projlinesvect = ones(1,length(cluster_dip_models));
  elseif strcmpi(options{2},'off') && strcmp(STUDY.etc.dipparams.centrline,'on')
      projlinesvect = zeros(1,length(cluster_dip_models));
      projlinesvect(centroidIdx) = 1;
  elseif strcmpi(options{2},'off')
      projlinesvect = zeros(1,length(cluster_dip_models));
  end
  options{2} = projlinesvect;

  dipplot(cluster_dip_models, options{:});
  set(fig_h, 'resize','on');

  %data cursor will show component label
  datacursormode on
  dcm = datacursormode(gcf);
  set(dcm,'Enable','on', 'UpdateFcn', @componentText)

end % multicolor. Supporting functions at end of file

% ========================================================================================

% STD_PLOTCOMPDIP - Commandline function, to visualizing cluster components dipoles.
%                   Displays the dipoles of specified cluster components with the cluster mean
%                   dipole on separate figures.
%                   To visualize dipoles they first must be stored in the EEG dataset structures
%                   using DIPFIT. Only components that have a dipole locations will be displayed,
%                   along with the cluster mean dipole in red.
% Usage:
%                   >> [STUDY] = std_plotcompdip(STUDY, ALLEEG, cluster, comps);
% Inputs:
%   STUDY      - EEGLAB STUDY set comprising some or all of the EEG datasets in ALLEEG.
%   ALLEEG     - global EEGLAB vector of EEG structures for the dataset(s) included in the STUDY.
%                     ALLEEG for a STUDY set is typically created using LOAD_ALLEEG.
%   cluster     - single cluster number.
%
% Optional inputs:
%   comps      - [numeric vector]  -> indices of the cluster components to plot.
%                       'all'                       -> plot all the components in the cluster {default: 'all'}.
%
% Outputs:
%   STUDY    - the input STUDY set structure modified with plotted cluster
%                     dipole mean, to allow quick replotting (unless cluster mean
%                     already existed in the STUDY).
%
%   Example:
%                         >> cluster = 4; comps= 1;
%                         >> [STUDY] = std_plotcompdip(STUDY,ALLEEG, cluster, comps);
%                    Plots component 1 dipole in blue with the cluster 4 mean dipole in red.
%
%  See also  pop_clustedit, dipfit, std_dipplot
%
% Authors:  Hilit Serby, Arnaud Delorme, Scott Makeig, SCCN, INC, UCSD, June, 2005

% Copyright (C) Hilit Serby, SCCN, INC, UCSD, June 08, 2005, hilit@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function STUDY = std_plotcompdip(STUDY, ALLEEG, cls, comp_ind, varargin)
if ~exist('cls')
    error('std_plotcompdip: you must provide a cluster number as an input.');
end
if isempty(cls)
   error('std_plotcompdip: you must provide a cluster number as an input.');
end
if nargin == 3 % no components indices were given
    % Default plot all components of the cluster
    [STUDY] = std_dipplot(STUDY, ALLEEG, 'clusters', cls);
    return
end

for ci = 1:length(comp_ind)
    abset = STUDY.datasetinfo(STUDY.cluster(cls).sets(1,comp_ind(ci))).index;
    comp = STUDY.cluster(cls).comps(comp_ind(ci));
    subject = STUDY.datasetinfo(STUDY.cluster(cls).sets(1,comp_ind(ci))).subject;
    if ~isfield(ALLEEG(abset), 'dipfit')
        warndlg2(['No dipole information available in dataset ' num2str(abset) ' , abort plotting'], 'Aborting plot dipoles');
        return;
    end
    if length(comp_ind) == 1 && isempty(ALLEEG(abset).dipfit.model(comp).posxyz)
        warndlg2(strvcat('There is no dipole information available in', ...
                       [ 'dataset ' num2str(abset) ' for this component, abort plotting']), 'Aborting plot dipoles');
        return;
    end
    if ~isfield(STUDY.cluster(cls),'dipole')
        STUDY = std_centroid(STUDY,ALLEEG, cls , 'dipole');
    elseif isempty(STUDY.cluster(cls).dipole)
        STUDY = std_centroid(STUDY,ALLEEG, cls , 'dipole');
    end
    comp_to_disp = [subject  ' / ' 'IC' num2str(comp) ];
    cluster_dip_models.posxyz = ALLEEG(abset).dipfit.model(comp).posxyz;
    cluster_dip_models.momxyz = ALLEEG(abset).dipfit.model(comp).momxyz;
    cluster_dip_models.rv = ALLEEG(abset).dipfit.model(comp).rv;
    cluster_dip_models(2) = STUDY.cluster(cls).dipole;
    if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
        if isfield(ALLEEG(abset).dipfit, 'hdmfile') %dipfit 2 spherical model
            load('-mat', ALLEEG(abset).dipfit.hdmfile);
            max_r = max(vol.r);
        else
            max_r = max(ALLEEG(abset).dipfit.vol.r);
        end
        dipplot(cluster_dip_models, 'sphere', max_r, 'mri', ALLEEG(abset).dipfit.mrifile,'coordformat', ALLEEG(abset).dipfit.coordformat , ...
           'normlen' ,'on', 'pointout' ,'on','color', {'b', 'r'}, 'dipnames', {comp_to_disp [ STUDY.cluster(cls).name ' mean' ] },...
            'spheres', 'on', 'verbose', 'off', varargin{:});
    else
       dipplot(cluster_dip_models, 'meshdata', ALLEEG(abset).dipfit.hdmfile, 'mri', ALLEEG(abset).dipfit.mrifile,'coordformat', ALLEEG(abset).dipfit.coordformat , ...
          'normlen' ,'on', 'pointout' ,'on','color', {'b', 'r'}, 'dipnames', {comp_to_disp [STUDY.cluster(cls).name ' mean']}, ...
          'spheres', 'off', 'verbose', 'off', varargin{:});
    end
    fig_h = gcf;
    set(fig_h,'Name', [subject ' / ' 'IC' num2str(comp) ', ' STUDY.cluster(cls).name],'NumberTitle','off');
end

% -----------------------
% load all dipoles and
% compute dipole centroid
% DEVELOPMENT: this function
% should be the only one to
% access dipole information
% -----------------------
function STUDY = std_centroid(STUDY,ALLEEG, clsind, tmp);

    for clust = 1:length(clsind)
        max_r = 0;
        len = length(STUDY.cluster(clsind(clust)).comps);
        tmppos = [ 0 0 0 ];
        tmpmom = [ 0 0 0 ];
        tmprv = 0;
        ndip = 0;
        for k = 1:len
            fprintf('.');
            comp  = STUDY.cluster(clsind(clust)).comps(k);
            abset = STUDY.cluster(clsind(clust)).sets(1,k);
            if ~isfield(ALLEEG(abset), 'dipfit')
               warndlg2(['No dipole information available in dataset ' num2str(abset) ], 'Aborting compute centroid dipole');
               return;
            end
            if ~isempty(ALLEEG(abset).dipfit.model(comp).posxyz)
                ndip   = ndip +1;
                posxyz = ALLEEG(abset).dipfit.model(comp).posxyz;
                momxyz = ALLEEG(abset).dipfit.model(comp).momxyz;
                if size(posxyz,1) == 2
                    if all(posxyz(2,:) == [ 0 0 0 ])
                        posxyz(2,:) = [];
                        momxyz(2,:) = [];
                    end
                end
                tmppos = tmppos + mean(posxyz,1);
                tmpmom = tmpmom + mean(momxyz,1);
                tmprv = tmprv + ALLEEG(abset).dipfit.model(comp).rv;
                if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
                   if isfield(ALLEEG(abset).dipfit, 'hdmfile') %dipfit 2 spherical model
                       load('-mat', ALLEEG(abset).dipfit.hdmfile);
                       max_r = max(max_r, max(vol.r));
                   else % old version of dipfit
                       max_r = max(max_r,max(ALLEEG(abset).dipfit.vol.r));
                   end
               end
            end
        end
        centroid{clust}.dipole.posxyz =  tmppos/ndip;
        centroid{clust}.dipole.momxyz =  tmpmom/ndip;
        centroid{clust}.dipole.rv =  tmprv/ndip;
        if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical') && (~isfield(ALLEEG(abset).dipfit, 'hdmfile')) %old dipfit
            centroid{clust}.dipole.maxr = max_r;
        end
        STUDY.cluster(clsind(clust)).dipole = centroid{clust}.dipole;
    end
    fprintf('\n');

% --------------------------------
% new function to compute centroid
% was programmed to debug the function
% above but is now used in the code
% --------------------------------
function dipole = computecentroid(alldipoles)

        max_r = 0;
        len = length(alldipoles);
        dipole.posxyz = [ 0 0 0 ];
        dipole.momxyz = [ 0 0 0 ];
        dipole.rv = 0;
        ndip = 0;
        count = 0;
        warningon = 1;
        for k = 1:len
            if size(alldipoles(k).posxyz,1) == 2
                if all(alldipoles(k).posxyz(2,:) == [ 0 0 0 ])
                    alldipoles(k).posxyz(2,:) = [];
                    alldipoles(k).momxyz(2,:) = [];
                end
            end
            if ~isempty(alldipoles(k).posxyz)
                dipole.posxyz = dipole.posxyz + mean(alldipoles(k).posxyz,1);
                dipole.momxyz = dipole.momxyz + mean(alldipoles(k).momxyz,1);
                dipole.rv     = dipole.rv     + alldipoles(k).rv;
                count = count+1;
            elseif warningon
                disp('Some components do not have dipole information');
                warningon = 0;
            end
        end
        dipole.posxyz = dipole.posxyz/count;
        dipole.momxyz = dipole.momxyz/count;
        dipole.rv     = dipole.rv/count;
        if isfield(alldipoles, 'maxr')
            dipole.maxr = alldipoles(1).max_r;
        end

function [cluster_dip_models, options] = dipgroups(ALLEEG, STUDY, cls, comp_to_disp, cluster_dip_models, options);

    % first, extract the subject number
    for n = 1:length(comp_to_disp)
        subjectnum(n,1) = str2num(comp_to_disp{n}(1:3));
    end

    % second, extract group info
    for n = 1:length(subjectnum)
        subj_group{n,1} = ALLEEG(1,subjectnum(n)).group;
    end

    % third, replace the group names with numbers
    for n = 1:length(subj_group)
        for m = 1:length(STUDY.group)
            if strcmp(subj_group{n,1}, STUDY.group{1,m})
                subj_groupnum(n,1) = m;
                break
            end
        end
    end

    % fourth, compute centroid for each group
    for n = 1:length(STUDY.group)
        samegroupIC = find(subj_groupnum==n);
        cluster_dip_models(1,length(subj_groupnum)+n) = computecentroid(cluster_dip_models(1, samegroupIC));
    end

    % fifth, use subj_groupnum as a type of dipole color

        %%%%%%%%%%%%%%%%%%%%% color list %%%%%%%%%%%%%%%%%%%%%
        % This color list was developed for std_envtopo
        % 16 colors names officially supported by W3C specification for HTML
        colors{1,1}  = [1 1 1];            % White
        colors{2,1}  = [1 1 0];            % Yellow
        colors{3,1}  = [1 0 1];            % Fuchsia
        colors{4,1}  = [1 0 0];            % Red
        colors{5,1}  = [0.75  0.75  0.75]; % Silver
        colors{6,1}  = [0.5 0.5 0.5];      % Gray
        colors{7,1}  = [0.5 0.5 0];        % Olive
        colors{8,1}  = [0.5 0 0.5];        % Purple
        colors{9,1}  = [0.5 0 0];          % Maroon
        colors{10,1} = [0 1 1];            % Aqua
        colors{11,1} = [0 1 0];            % Lime
        colors{12,1} = [0 0.5 0.5];        % Teal
        colors{13,1} = [0 0.5 0];          % Green
        colors{14,1} = [0 0 1];            % Blue
        colors{15,1} = [0 0 0.5];          % Navy
        colors{16,1} = [0 0 0];            % Black
        % Silver is twice brighter because silver is used for a background color
        colors{5,1} = [0.875 0.875 0.875];
        % Choosing and sorting 12 colors for line plot, namely Red, Blue, Green, Fuchsia, Lime, Aqua, Maroon, Olive, Purple, Teal, Navy, and Gray
        selectedcolors = colors([4 13 14 3 11 10 9 7 8 12 15 6]);

    % determine the new dip colors
    for n = 1:length(subj_groupnum)
        dip_color{1,n}=selectedcolors{subj_groupnum(n,1)+1};
    end
    for n = 1:length(STUDY.group)
        dip_color{1,end+1}= selectedcolors{n+1};
    end

    for n = 1:length(options)
        if      strcmp(options{1,n}, 'color')
            options{1,n+1} = dip_color;
        elseif  strcmp(options{1,n}, 'dipnames')
            dipnames = options{1,n+1};
            for m = 1:length(STUDY.group)
                dipnames{1,length(subj_groupnum)+m}= [STUDY.group{1,m} ' mean'];
            end
            options{1,n+1} = dipnames;
        end
    end

% ========================================================================================
% multicolor mode support functions
function str = componentText(~,obj)
% look up component name--used in data cursor callback in 'multicolor' mode
try
  str = obj.Target.UserData.name;
  h=findall(gcf,'type','hggroup');
  h(1).FontSize = 14;
catch
  str = '';
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [out] = evap_22(p1,p2,S,Ep,dt)
%evap_22 3-part piece-wise evaporation

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Threshold-based evaporation rate
% Constraints:  f <= S/dt
% @(Inputs):    p1   - wilting point [mm]
%               p2   - 2nd (lower) threshold [mm]
%               S    - current storage [mm]
%               Ep   - potential evapotranspiration rate [mm/d]
%               dt   - time step size [d]

out = min(S/dt,max(0,min((S-p1)./(p2-p1).*Ep,Ep)));

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function ulimit = SlopeLimitLin(ul,xl,vm1,v0,vp1);

% function ulimit = SlopeLimitLin(ul,xl,vm1,v0,vp1);
% Purpose: Apply slopelimited on linear function ul(Np,1) on x(Np,1)
%          (vm1,v0,vp1) are cell averages left, center, and right

Globals1D;

% Compute various geometric measures
ulimit = ul; h = xl(Np,:)-xl(1,:);
x0 = ones(Np,1)*(xl(1,:) + h/2);

hN = ones(Np,1)*h;

% Limit function
ux = (2./hN).*(Dr*ul);

ulimit = ones(Np,1)*v0+(xl-x0).*(ones(Np,1)*minmod([ux(1,:);(vp1-v0)./h;(v0-vm1)./h]));
%ulimit = ones(Np,1)*v0+(xl-x0).*(ones(Np,1)*minmodB([ux(1,:);(vp1-v0)./h;(v0-vm1)./h],1e-2,h));
return

repo: Aero-matlab


text:
function pn = line_par_point_near_3d ( f, g, h, x0, y0, z0, p )

%*****************************************************************************80
%
%% LINE_PAR_POINT_NEAR_3D: nearest point on parametric line to given point, 3D.
%
%  Discussion:
%
%    The parametric form of a line in 3D is:
%
%      X = X0 + F * T
%      Y = Y0 + G * T
%      Z = Z0 + H * T
%
%    We may normalize by choosing F*F + G*G + H*H = 1, and F nonnegative.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer, John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983,
%    ISBN: 0408012420.
%
%  Parameters:
%
%    Input, real F, G, H, X0, Y0, Z0, the parametric
%    line parameters.
%
%    Input, real P(3,1), the point whose distance from the line is
%    to be measured.
%
%    Output, real PN(3,1), the point on the parametric line which
%    is nearest to P.
%
  t = ( f * ( p(1) - x0 ) + g * ( p(2) - y0 ) + h * ( p(3) - z0 ) ) ...
    / ( f * f + g * g + h * h );

  pn(1,1) = x0 + t * f;
  pn(2,1) = y0 + t * g;
  pn(3,1) = z0 + t * h;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Faces, iFacesRemove] = tess_threshold(Vertices, Faces, threshArea, threshRatio, threshAngle, threshEdge)
% TESS_THRESHOLD: Detect pathological triangles in a surface mesh.
%
% INPUTS:
%    - Vertices    : [Nvert x 3] surface vertices
%    - Faces       : [Nfaces x 3] surface triangles
%    - threshArea  : Detect large triangles (area > thresh * std)
%    - threshRatio : Detect asymetric triangles (ratio perimeter/area > thresh * std)
%    - threshAngle : Detect triangles with angles that are too open (angle in degrees > thresh)

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

% ===== PARSE INPUTS =====
if (nargin < 6) || isempty(threshEdge)
    threshEdge = [];
end
if (nargin < 5) || isempty(threshAngle)
    threshAngle = [];
end
if (nargin < 4) || isempty(threshRatio)
    threshRatio = [];
end
if (nargin < 3)
    threshArea = [];
end
iFacesRemoveArea  = [];
iFacesRemoveRatio = [];
iFacesRemoveAngle = [];
iFacesRemoveEdge  = [];

% ===== COMPUTE SURFACE STATISTICS =====
% Triangles area
triArea = tess_area(Vertices, Faces);
% Compute the vector of each edge
v1 = Vertices(Faces(:,1),:) - Vertices(Faces(:,2),:);
v2 = Vertices(Faces(:,1),:) - Vertices(Faces(:,3),:);
v3 = Vertices(Faces(:,2),:) - Vertices(Faces(:,3),:);

% ===== THRESHOLD: AREA =====
% Detect the faces that have an area above the threshold
if ~isempty(threshArea) && (threshArea > 0)
    iFacesRemoveArea = find(triArea - mean(triArea) > threshArea * std(triArea));
end

% ===== THRESHOLD: PERIMETER/AREA =====
if ~isempty(threshRatio) && (threshRatio > 0)
    % Compute perimeter again
    triPerimeter = tess_perimeter(Vertices, Faces);
    % Ratio perimeter / area
    ratio = (triPerimeter ./ triArea);
    % Detect the Faces that have an area above the threshold
    iFacesRemoveRatio = find(ratio - mean(ratio) > threshRatio * std(ratio));
end

% ===== THRESHOLD: ANGLE =====
if ~isempty(threshAngle) && (threshAngle > 0)
    % Compute the angle between all the vectors
    maxAngle = zeros(size(Vertices,1),1);
    for i = 1:size(v1,1)
        maxAngle(i) = max([atan2(norm(cross(v1(i,:),v2(i,:))), dot(v1(i,:),v2(i,:))), ...
                        atan2(norm(cross(v1(i,:),v3(i,:))), dot(v1(i,:),v3(i,:))), ...
                        atan2(norm(cross(v2(i,:),v3(i,:))), dot(v2(i,:),v3(i,:)))]);
    end
    % Convert to degrees
    maxAngle = maxAngle / 2 / pi * 360;
    % Detect the Faces that have an area above the threshold
    iFacesRemoveAngle = find(maxAngle > threshAngle);
end

% ===== THRESHOLD: EDGE LENGTH =====
if ~isempty(threshEdge) && (threshEdge > 0)
    % Compute the length of all the edges
    edgeLength = sqrt(v1.^2 + v2.^2 + v3.^2);
    % Split long edges
    iFacesRemoveEdge = find(edgeLength - mean(edgeLength) > threshEdge * std(edgeLength));
end

% List of faces to remove
iFacesRemove = [iFacesRemoveArea(:); iFacesRemoveRatio(:); iFacesRemoveAngle(:)];
% Keep only the good faces
if ~isempty(iFacesRemove)
     Faces(iFacesRemove,:) = [];
end

author: brainstorm-tools


text:
function xyb = Xb(s)

x = s ;

if s <= 0.5
    y = -s ;
elseif s>0.5
    y = s-1 ;
end

xyb = [x ; y] ;

author: Sable

repo: mcbench-benchmarks


text:
%SETCOST Reset classification cost matrix of mapping
%
%   W = SETCOST(W,COST,LABLIST)
%
% The classification cost matrix of the dataset W is reset to COST.
% W has to be a trained classifier. COST should have size [C,C+1],
% if C is the number of classes assigned by W.
% COST(I,J) are the costs of classifying an object of class I
% as class J. Column C+1 generates an alternative reject class and
% may be omitted, yielding a size of [C,C].
% An empty cost matrix, COST = [] (default) is interpreted as
% COST = ONES(C) - EYE(C) (identical costs of misclassification).
%
% In LABLIST the corresponding class labels may be supplied.
% LABLIST may have only class names of the existing labels assigned
% by W, stored in W.LABELS.

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [dims,scales,bpp,endian,datatype] = read_avw_hdr(fname)
% [dims,scales,bpp,endian] = READ_AVW_HDR(fname)
%
%  Extracts the 4 dimensions (dims),
%   4 scales (scales) and bytes per pixel (bpp) for voxels
%   contained in the Analyse header file (fname)
%   Also returns endian = 'l' for little-endian or 'b' for big-endian
%
%  See also: READ_AVW, READ_AVW_IMG, SAVE_AVW, SAVE_AVW_HDR, SAVE_AVW_IMG

% remove extension if it exists
if ( (length(findstr(fname,'.hdr'))>0) | ...
        (length(findstr(fname,'.img')>0)) ),
  fname=fname(1:(length(fname)-4));
end
fnhdr=strcat(fname,'.hdr');

% open file in big-endian
endian='b';
fid=fopen(fnhdr,'r','b');
testval = fread(fid,1,'int32');
% check if this gives the correct header size - if not use little-endian
if (testval~=348),
  fclose(fid);
  fid=fopen(fnhdr,'r','l');
  endian='l';
  testval = fread(fid,1,'int32');
  if (testval~=348),
    disp('Can not read this file format');
    return;
  end
end
        % ditch the remaining initial header stuff
  dummy=fread(fid,36,'char');
        % ditch dim[0] = No. dimensions
  dummy=fread(fid,1,'int16');
  dims=fread(fid,4,'int16');
  dummy=fread(fid,3,'int16');
  dummy=fread(fid,14,'char');
  datatype=fread(fid,1,'int16');
  bpp=fread(fid,1,'int16');
  dummy=fread(fid,2,'char');
  dummy=fread(fid,1,'float');
  scales=fread(fid,4,'float');
fclose(fid);
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function realtest_AD1()
% Test AD for a real optimization problem on a product manifold (struct)

    % Verify that Manopt was indeed added to the Matlab path.
    if isempty(which('spherefactory'))
        error(['You should first add Manopt to the Matlab path.\n' ...
		       'Please run importmanopt.']);
    end

    % Verify that the deep learning tool box was installed
    assert(exist('dlarray', 'file') == 2, ['Deep learning tool box is '...
    'needed for automatic differentiation.\n Please install the'...
    'latest version of the deep learning tool box and \nupgrade to Matlab'...
    ' R2021b if possible.'])

    % Generate the problem data.
    n = 100;
    A = randn(n);
    A = .5*(A+A');

    % Create the product manifold
    S = spherefactory(n);
    manifold.x = S;
    manifold.y = S;
    problem.M = productmanifold(manifold);

    % Define the problem cost function
    problem.cost  = @(X) -X.x'*(A*X.y);

    % Define the gradient and the hessian via automatic differentiation
    problem = manoptAD(problem);

    % Numerically check gradient and Hessian consistency.
    figure;
    checkgradient(problem);
    figure;
    checkhessian(problem);

    % Solve.
    [x, xcost, info] = trustregions(problem);          %#ok<ASGLU>

    % Test
    ground_truth = svd(A);
    distance = abs(ground_truth(1) - (-problem.cost(x)));
    fprintf('The distance between the ground truth and the solution is %e \n',distance);

end

author: NicolasBoumal

repo: manopt


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%  	System response to sinusoidal inputs

% which element of vector a has the closest value to z
a=-4:4
z=2.3;
a-z
abs(a-z)
[m,i]=min(abs(a-z))
a(i)

%  	System response
num=[3 4 2];
den=[1 1 3];
[H,w]=freqs(num,den);
w0=1;
[m,i]=min(abs(w-w0));
w(i)
Hw0=H(i)
mag=abs(Hw0)
phas=angle(Hw0)
t=0:.1:30;
y1=3*mag*cos(w(i)*t+pi/4+phas)
plot(t,y1)
legend('Output  y(t)')
ylim([-6 8]);


text:
function boolVal=isImag(x)
%%ISIMAG Returns true if x is a purely imaginary vector or matrix. That
%        is, it is complex and all real elements are exactly zero.
%
%INPUT: x A scalar or matrix. This must be of a numeric type; that is, not
%         a cell array, a string, etc.
%
%OUTPUTS: boolVal This is true is x is a complex matrix (isreal is false)
%                 and all real components are zero.
%
%December 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(isreal(x))
   boolVal=false;
   return;
end

if(all(real(x(:))==0))
    boolVal=true;
else
    boolVal=false;
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%ELLIPSE2POLY  Approximates an elliptic arc with a polyline
%
%     pts = cv.ellipse2Poly(center, axes)
%     [...] = cv.ellipse2Poly(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __center__ Center of the arc `[x,y]`.
% * __axes__ Half of the size of the ellipse main axes `[a,b]`. See cv.ellipse
%   for details.
%
% ## Output
% * __pts__ Output vector of polyline vertices. An Nx2 numeric matrix
%   `[x y; ...]`.
%
% ## Options
% * __Angle__ Rotation angle of the ellipse in degrees. See cv.ellipse for
%   details. default 0.
% * __StartAngle__ Starting angle of the elliptic arc in degrees. default 0
% * __EndAngle__ Ending angle of the elliptic arc in degrees. default 360
% * __Delta__ Angle between the subsequent polyline vertices. It defines the
%   approximation accuracy. default 5.
%
% The function cv.ellipse2Poly computes the vertices of a polyline that
% approximates the specified elliptic arc. It is used by cv.ellipse. If
% `StartAngle` is greater than `EndAngle`, they are swapped.
%
% See also: cv.ellipse
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%  Save NIFTI or ANALYZE dataset that is loaded by "load_untouch_nii.m".
%  The output image format and file extension will be the same as the
%  input one (NIFTI.nii, NIFTI.img or ANALYZE.img). Therefore, any file
%  extension that you specified will be ignored.
%
%  Usage: save_untouch_nii(nii, filename)
%
%  nii - nii structure that is loaded by "load_untouch_nii.m"
%
%  filename  - 	NIFTI or ANALYZE file name.
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function save_untouch_nii(nii, filename)

   if ~exist('nii','var') | isempty(nii) | ~isfield(nii,'hdr') | ...
	~isfield(nii,'img') | ~exist('filename','var') | isempty(filename)

      error('Usage: save_untouch_nii(nii, filename)');
   end

   if ~isfield(nii,'untouch') | nii.untouch == 0
      error('Usage: please use ''save_nii.m'' for the modified structure.');
   end

   if isfield(nii.hdr.hist,'magic') & strcmp(nii.hdr.hist.magic(1:3),'ni1')
      filetype = 1;
   elseif isfield(nii.hdr.hist,'magic') & strcmp(nii.hdr.hist.magic(1:3),'n+1')
      filetype = 2;
   else
      filetype = 0;
   end

   [p,f] = fileparts(filename);
   fileprefix = fullfile(p, f);

   write_nii(nii, filetype, fileprefix);

%   %  So earlier versions of SPM can also open it with correct originator
 %  %
  % if filetype == 0
   %   M=[[diag(nii.hdr.dime.pixdim(2:4)) -[nii.hdr.hist.originator(1:3).*nii.hdr.dime.pixdim(2:4)]'];[0 0 0 1]];
    %  save(fileprefix, 'M');
%   elseif filetype == 1
 %     M=[];
  %    save(fileprefix, 'M');
   %end

   return					% save_untouch_nii

%-----------------------------------------------------------------------------------
function write_nii(nii, filetype, fileprefix)

   hdr = nii.hdr;

   if isfield(nii,'ext') & ~isempty(nii.ext)
      ext = nii.ext;
      [ext, esize_total] = verify_nii_ext(ext);
   else
      ext = [];
   end

   switch double(hdr.dime.datatype),
   case   1,
      hdr.dime.bitpix = int16(1 ); precision = 'ubit1';
   case   2,
      hdr.dime.bitpix = int16(8 ); precision = 'uint8';
   case   4,
      hdr.dime.bitpix = int16(16); precision = 'int16';
   case   8,
      hdr.dime.bitpix = int16(32); precision = 'int32';
   case  16,
      hdr.dime.bitpix = int16(32); precision = 'float32';
   case  32,
      hdr.dime.bitpix = int16(64); precision = 'float32';
   case  64,
      hdr.dime.bitpix = int16(64); precision = 'float64';
   case 128,
      hdr.dime.bitpix = int16(24); precision = 'uint8';
   case 256
      hdr.dime.bitpix = int16(8 ); precision = 'int8';
   case 512
      hdr.dime.bitpix = int16(16); precision = 'uint16';
   case 768
      hdr.dime.bitpix = int16(32); precision = 'uint32';
   case 1024
      hdr.dime.bitpix = int16(64); precision = 'int64';
   case 1280
      hdr.dime.bitpix = int16(64); precision = 'uint64';
   case 1792,
      hdr.dime.bitpix = int16(128); precision = 'float64';
   otherwise
      error('This datatype is not supported');
   end

%   hdr.dime.glmax = round(double(max(nii.img(:))));
 %  hdr.dime.glmin = round(double(min(nii.img(:))));

   if filetype == 2
      fid = fopen(sprintf('%s.nii',fileprefix),'w');

      if fid < 0,
         msg = sprintf('Cannot open file %s.nii.',fileprefix);
         error(msg);
      end

      hdr.dime.vox_offset = 352;

      if ~isempty(ext)
         hdr.dime.vox_offset = hdr.dime.vox_offset + esize_total;
      end

      hdr.hist.magic = 'n+1';
      save_untouch_nii_hdr(hdr, fid);

      if ~isempty(ext)
         save_nii_ext(ext, fid);
      end
   elseif filetype == 1
      fid = fopen(sprintf('%s.hdr',fileprefix),'w');

      if fid < 0,
         msg = sprintf('Cannot open file %s.hdr.',fileprefix);
         error(msg);
      end

      hdr.dime.vox_offset = 0;
      hdr.hist.magic = 'ni1';
      save_untouch_nii_hdr(hdr, fid);

      if ~isempty(ext)
         save_nii_ext(ext, fid);
      end

      fclose(fid);
      fid = fopen(sprintf('%s.img',fileprefix),'w');
   else
      fid = fopen(sprintf('%s.hdr',fileprefix),'w');

      if fid < 0,
         msg = sprintf('Cannot open file %s.hdr.',fileprefix);
         error(msg);
      end

      save_untouch0_nii_hdr(hdr, fid);

      fclose(fid);
      fid = fopen(sprintf('%s.img',fileprefix),'w');
   end

   ScanDim = double(hdr.dime.dim(5));		% t
   SliceDim = double(hdr.dime.dim(4));		% z
   RowDim   = double(hdr.dime.dim(3));		% y
   PixelDim = double(hdr.dime.dim(2));		% x
   SliceSz  = double(hdr.dime.pixdim(4));
   RowSz    = double(hdr.dime.pixdim(3));
   PixelSz  = double(hdr.dime.pixdim(2));

   x = 1:PixelDim;

   if filetype == 2 & isempty(ext)
      skip_bytes = double(hdr.dime.vox_offset) - 348;
   else
      skip_bytes = 0;
   end

   if double(hdr.dime.datatype) == 128

      %  RGB planes are expected to be in the 4th dimension of nii.img
      %
      if(size(nii.img,4)~=3)
         error(['The NII structure does not appear to have 3 RGB color planes in the 4th dimension']);
      end

      nii.img = permute(nii.img, [4 1 2 3 5 6 7 8]);
   end

   %  For complex float32 or complex float64, voxel values
   %  include [real, imag]
   %
   if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
      real_img = real(nii.img(:))';
      nii.img = imag(nii.img(:))';
      nii.img = [real_img; nii.img];
   end

   if skip_bytes
      fwrite(fid, zeros(1,skip_bytes), 'uint8');
   end

   fwrite(fid, nii.img, precision);
%   fwrite(fid, nii.img, precision, skip_bytes);        % error using skip
   fclose(fid);

   return;					% write_nii


text:
%  Verify NIFTI header extension to make sure that each extension section
%  must be an integer multiple of 16 byte long that includes the first 8
%  bytes of esize and ecode. If the length of extension section is not the
%  above mentioned case, edata should be padded with all 0.
%
%  Usage: [ext, esize_total] = verify_nii_ext(ext)
%
%  ext - Structure of NIFTI header extension, which includes num_ext,
%       and all the extended header sections in the header extension.
%       Each extended header section will have its esize, ecode, and
%       edata, where edata can be plain text, xml, or any raw data
%       that was saved in the extended header section.
%
%  esize_total - Sum of all esize variable in all header sections.
%
%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function [ext, esize_total] = verify_nii_ext(ext)

   if ~isfield(ext, 'section')
      error('Incorrect NIFTI header extension structure.');
   elseif ~isfield(ext, 'num_ext')
      ext.num_ext = length(ext.section);
   elseif ~isfield(ext, 'extension')
      ext.extension = [1 0 0 0];
   end

   esize_total = 0;

   for i=1:ext.num_ext
      if ~isfield(ext.section(i), 'ecode') | ~isfield(ext.section(i), 'edata')
         error('Incorrect NIFTI header extension structure.');
      end

      ext.section(i).esize = ceil((length(ext.section(i).edata)+8)/16)*16;
      ext.section(i).edata = ...
	[ext.section(i).edata ...
	 zeros(1,ext.section(i).esize-length(ext.section(i).edata)-8)];
      esize_total = esize_total + ext.section(i).esize;
   end

   return                                       % verify_nii_ext

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function outstr = mets2str(str)
% Creates strings for message boxes
%
% USAGE:
%
%    outstr = mets2str(str)
%
% INPUT:
%    str:       struct of names
%
% OUTPUT:
%    outstr:    'act', 'acon-C', 'pep'
%
% .. Author: - Stefan G. Thorleifsson 2011
%
% .. rBioNet is published under GNU GENERAL PUBLIC LICENSE 3.0+
% .. Thorleifsson, S. G., Thiele, I., rBioNet: A COBRA toolbox extension for
% .. reconstructing high-quality biochemical networks, Bioinformatics, Accepted.
% .. rbionet@systemsbiology.is

S = size(str);
if S(1) <= S(2)
    m = S(2);
else
    m = S(1);
end

outstr = [];

for i = 1:m
    if i == 1
        outstr = str{i};
    else
        outstr = [outstr ', ' str{i}];
    end
end

author: opencobra

repo: cobratoolbox


text:
function S = L0Restoration(Im, kernel, lambda, kappa)
%%
% Image restoration with L0 prior
% The objective function:
% S^* = argmin ||I*k - B||^2 + lambda |\nabla I|_0
%% Input:
% @Im: Blurred image
% @kernel: blur kernel
% @lambda: weight for the L0 prior
% @kappa: Update ratio in the ADM
%% Output:
% @S: Latent image
%
% The Code is created based on the method described in the following paper
%   [1] Jinshan Pan, Zhe Hu, Zhixun Su, and Ming-Hsuan Yang,
%        Deblurring Text Images via L0-Regularized Intensity and Gradient
%        Prior, CVPR, 2014.
%   [2] Li Xu, Cewu Lu, Yi Xu, and Jiaya Jia. Image smoothing via l0 gradient minimization.
%        ACM Trans. Graph., 30(6):174, 2011.
%
%   Author: Jinshan Pan (sdluran@gmail.com)
%   Date  : 05/18/2014

if ~exist('kappa','var')
    kappa = 2.0;
end
%% pad image
H = size(Im,1);    W = size(Im,2);
Im = wrap_boundary_liu(Im, opt_fft_size([H W]+size(kernel)-1));
%%
S = Im;
betamax = 1e5;
fx = [1, -1];
fy = [1; -1];
[N,M,D] = size(Im);
sizeI2D = [N,M];
otfFx = psf2otf(fx,sizeI2D);
otfFy = psf2otf(fy,sizeI2D);
%%
KER = psf2otf(kernel,sizeI2D);
Den_KER = abs(KER).^2;
%%
Denormin2 = abs(otfFx).^2 + abs(otfFy ).^2;
if D>1
    Denormin2 = repmat(Denormin2,[1,1,D]);
    KER = repmat(KER,[1,1,D]);
    Den_KER = repmat(Den_KER,[1,1,D]);
end
Normin1 = conj(KER).*fft2(S);
%%
beta = 2*lambda;
while beta < betamax
    Denormin   = Den_KER + beta*Denormin2;
    h = [diff(S,1,2), S(:,1,:) - S(:,end,:)];
    v = [diff(S,1,1); S(1,:,:) - S(end,:,:)];
    if D==1
        t = (h.^2+v.^2)<lambda/beta;
    else
        t = sum((h.^2+v.^2),3)<lambda/beta;
        t = repmat(t,[1,1,D]);
    end
    h(t)=0; v(t)=0;
    Normin2 = [h(:,end,:) - h(:, 1,:), -diff(h,1,2)];
    Normin2 = Normin2 + [v(end,:,:) - v(1, :,:); -diff(v,1,1)];
    FS = (Normin1 + beta*fft2(Normin2))./Denormin;
    S = real(ifft2(FS));
    beta = beta*kappa;
end
S = S(1:H, 1:W, :);
end

author: IVRL

repo: Kernel-Modeling-Super-Resolution


text:
function [ w, xy ] = triangle_unit_o03 ( )

%*****************************************************************************80
%
%% TRIANGLE_UNIT_O03 returns a 3 point quadrature rule for the unit triangle.
%
%  Discussion:
%
%    This rule is precise for monomials through degree 2.
%
%    The integration region is:
%
%      0 <= X
%      0 <= Y
%      X + Y <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carlos Felippa,
%    A compendium of FEM integration formulas for symbolic work,
%    Engineering Computation,
%    Volume 21, Number 8, 2004, pages 867-890.
%
%  Parameters:
%
%    Output, real W(3), the weights.
%
%    Output, real XY(2,3), the abscissas.
%
  w(1:3,1) = [ ...
    0.33333333333333333333, ...
    0.33333333333333333333, ...
    0.33333333333333333333 ];

  xy(1:2,1:3) = [ ...
    0.66666666666666666667,  0.16666666666666666667; ...
    0.16666666666666666667,  0.66666666666666666667; ...
    0.16666666666666666667,  0.16666666666666666667 ]';
  return
end


text:
function analyse_eval_lin_fusion(key_dev,ndx_eval,dev_scores_obj_array,eval_scores_obj_array,dev_scr_lin,eval_scr_lin)
% Does some analyses on evaluation scores for all systems and all
% fusions.  For each system, it calculates the estimated percentage
% of targets in the database.  It also calculates a KL divergence
% matrix of all systems against all systems.
% Inputs:
%   key_dev: The Key for the dev trials.
%   ndx_eval: The Ndx describing the eval trials.
%   dev_scores_obj_array: An array of Scores objects containing
%     system scores for the dev trials.
%   eval_scores_obj_array: An array of Scores objects containing
%     system scores for the eval trials.
%   dev_scr_lin: The dev scores resulting from the linear fusion.
%   eval_scr_lin: The eval scores resulting from the linear fusion.

assert(nargin==6)
assert(isa(key_dev,'Key'))
assert(isa(ndx_eval,'Ndx'))
assert(isa(dev_scr_lin,'Scores'))
assert(isa(eval_scr_lin,'Scores'))
assert(key_dev.validate())
assert(ndx_eval.validate())
assert(dev_scr_lin.validate())
assert(eval_scr_lin.validate())

eval_stacked_scores = stackScores(ndx_eval,eval_scores_obj_array);
trialmask = ndx_eval.trialmask;
eval_lin_fusion = eval_scr_lin.scoremat(trialmask(:))';
eval_scores = [eval_stacked_scores;eval_lin_fusion];

ptar = estimate_target_proportions(eval_scores,0.08);
logprint(Logger.Info,'\nSystem Ptar estimates:\n');
logprint(Logger.Info,'%s',sprintfmatrix(100*ptar'));

ptar_best = ptar(end);

dev_stacked_scores = stackScores(key_dev,dev_scores_obj_array);
trialmask = key_dev.to_ndx().trialmask;
dev_lin_fusion = dev_scr_lin.scoremat(trialmask(:))';
dev_scores = [dev_stacked_scores;dev_lin_fusion];

devKL = 100*kldivergence_of_systems(dev_scores,ptar_best,key_dev);
evalKL = 100*kldivergence_of_systems(eval_scores,ptar_best);
logprint(Logger.Info,'\nDev score divergence:\n%s',sprintfmatrix(devKL));
logprint(Logger.Info,'\nEval score divergence:\n%s',sprintfmatrix(evalKL));
end

author: nesl


text:
function exact = p20_exact ( )

%*****************************************************************************80
%
%% P20_EXACT returns the estimated integral for problem 20.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the estimated value of the integral.
%
  exact = 1.5643964440690497731;

  return
end


text:
function [ids, ids_com] = pickRandomSubsetIndex(n, k)
	%% ================== File info ==========================
	% Author		: Tiep Vu (http://www.personal.psu.edu/thv102/)
	% Time created	: Tue Jan 26 22:59:40 2016
	% Last modified	: Tue Jan 26 22:59:41 2016
	% Description	: pick a k-element subset of the set 1: n
	% 	INPUT:
	%		n: number of elements
	%		k: number of picked elements
	% 	OUTPUT:
	%		ids: vector of picked indices
	% 		ids_com: vector of unpicked indices
	%
	%% ================== end File info ==========================
    id_mixed = randperm(n);
    ids = id_mixed(1: k);
    if nargout == 2
    	ids_com = setdiff(1:n, ids);
    end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function smooth = score_smooth(score)
score = score';
tmp1 = [score(2:end),score(end)];
tmp2 = [score(1), score(1:end-1)];
smooth = (score +  tmp1 + 0*tmp2)/2;
end

author: wanglimin

repo: UntrimmedNet


text:
function whitehKernDisplay(kern, spacing)

% WHITEHKERNDISPLAY Display parameters of the WHITEH kernel.
% FORMAT
% DESC displays the parameters of the white noise
% kernel and the kernel type to the console.
% ARG kern : the kernel to display.
%
% FORMAT does the same as above, but indents the display according
% to the amount specified.
% ARG kern : the kernel to display.
% ARG spacing : how many spaces to indent the display of the kernel by.
%
% SEEALSO : whiteKernParamInit, modelDisplay, kernDisplay
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

if nargin > 1
  spacing = repmat(32, 1, spacing);
else
  spacing = [];
end
spacing = char(spacing);
fprintf(spacing);
fprintf('White Noise Variance: %2.4f\n', kern.variance)

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function value = i4_max ( a, b )

%*****************************************************************************80
%
%% I4_MAX returns the maximum of two I4's.
%
%  Discussion:
%
%    An I4 is an integer value.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer A, B, values to compare.
%
%    Output, real VALUE, the maximum of A and B.
%
  a = floor ( a );
  b = floor ( b );

  if ( a < b )
    value = b;
  else
    value = a;
  end

  return
end


text:
% This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% Copyright: ETH Zurich, 2014
% Author: Romano Cicchetti

function [ FSMs ] = generateFSMs(clusters, numOfClusters, clusterSizes, numOfFSMs, fitnessWeights, maxNumOfStates)

    % generate finite state machines that correspond to potential appliances
    % select the numOfFSMs finite state machines with the highest fitness
    % values

    population = [];
    for i = 2:maxNumOfStates
        population = [population; combnk(1:numOfClusters, i)];
    end
    fitnessValues = computeFitness(clusters, clusterSizes, population, numOfClusters, fitnessWeights);
    [~, sortedIdx] = sort(fitnessValues);
    FSMs = population(sortedIdx(1:numOfFSMs), :);

end

author: beckel

repo: nilm-eval

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


text:
function [m_CO2, y_CO2] = DuanSun(T, p, m_salt)
% This function solves the Duan-Sun model to calculate the solubility of
% CO2 in water and in the brine and the mole fraction of CO2 in the gas phase
% T is the temperature in [K]
% p is the pressure in [Pa]
% m is the salt molality in [mol/kg]
% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file
% R = 0.08314467; % [bar.L/(mol.K)]
% Written by Ali A. Eftekhari at TU Delft

P = p/1e5; % convert Pa to bar
% Constants (Table 2):
% constants for the calculation of chemical potential of CO2 in the liquid phase
c_mu_CO2_L = [28.9447706, -0.0354581768, -4770.67077, 1.02782768e-5, ...
		  33.8126098, 9.04037140e-3, -1.14934031e-3, -0.307405726, ...
		  -0.0907301486, 9.32713393e-4, 0.0];
% Constants for the calculation of second and third order interaction parameters
c_lambda_CO2_Na = [-0.411370585, 6.07632013e-4, 97.5347708, 0.0, 0.0, ...
				   0.0, 0.0, -0.0237622469, 0.0170656236, 0.0, 1.41335834e-5];
c_zeta_CO2_NaCl = [3.36389723e-4, -1.98298980e-5, 0.0, 0.0, 0.0, ...
				    0.0, 0.0, 2.12220830e-3, -5.24873303e-3, 0.0, 0.0];

% using the above constants (Eq. 7)
par = @(c)(c(1)+c(2)*T+c(3)/T+c(4)*T^2+c(5)/(630-T)+c(6)*P+ ...
		   c(7)*P*log(T)+c(8)*P/T+c(9)*P/(630-T)+ ...
		   c(10)*P^2/(630-T)^2+c(11)*T*log(P));

mu_CO2_L = par(c_mu_CO2_L);
lambda_CO2_Na = par(c_lambda_CO2_Na);
zeta_CO2_NaCl = par(c_zeta_CO2_NaCl);

% activity coefficient of CO2 (Eq. 5)
m_Na = m_salt;
m_Cl = m_salt;
log_gama_CO2 = 2*lambda_CO2_Na*m_Na+zeta_CO2_NaCl*m_Na*m_Cl;

% calculation of the fugacity of CO2 using the equation of Duan
a = [8.99288497e-2, -4.94783127e-1, 4.77922245e-2, 1.03808883e-2, ...
     -2.82516861e-2, 9.49887563e-2, 5.20600880e-4, -2.93540971e-4, ...
     -1.77265112e-3, -2.51101973e-5, 8.93353441e-5, 7.88998563e-5, ...
     -1.66727022e-2, 1.39800000,     2.96000000e-2];
% CO2 critical properties
Tc = 304.25; % [K]
Pc = 7.39e6; % [Pa]
Vc = 8.314467*Tc/Pc;
Tr = T/Tc;
Pr = P*1.0e5/Pc;
Z = @(Vr)(1+(a(1)+a(2)/Tr^2+a(3)/Tr^3)/Vr + ...
          (a(4)+a(5)/Tr^2+a(6)/Tr^3)/Vr^2 + (a(7)+a(8)/Tr^2+a(9)/Tr^3)/Vr^4 ...
          + (a(10)+a(11)/Tr^2+a(12)/Tr^3)/Vr^5 + ...
          a(13)/(Tr^3*Vr^2)*(a(14)+a(15)/Vr^2)*exp(-a(15)/Vr^2));
Vr_guess = (8.314*T*0.9/(P*1.0e5))/Vc;
F = @(Vr)(Z(Vr)/Vr-Pr/Tr);
Vr = fzero(F, Vr_guess);

% fugacity coefficient of pure CO2
Z_num = Z(Vr);
log_phi_CO2 = Z_num - 1.0 - log(Z_num) + (a(1)+a(2)/Tr^2+a(3)/Tr^3)/Vr + ...
			(a(4)+a(5)/Tr^2+a(6)/Tr^3)/(2.0*Vr^2) + ...
			(a(7)+a(8)/Tr^2+a(9)/Tr^3)/(4.0*Vr^4) + ...
			(a(10)+a(11)/Tr^2+a(12)/Tr^3)/(5.0*Vr^5) + ...
			a(13)/(2.0*Tr^3*a(15))*(a(14)+1.0-(a(14)+1.0+a(15)/Vr^2)*exp(-a(15)/Vr^2));

% pure water vapor pressure
c = [-38.640844, 5.8948420, 59.876516, 26.654627, 10.637097];
Pcw = 220.85; % [bar]
Tcw = 647.29; % [K]
t = (T-Tcw)/Tcw;
P_water = Pcw*T/Tcw*(1+c(1)*(-t)^1.9+c(2)*t+c(3)*t^2+c(4)*t^3+c(5)*t^4);
y_CO2 = (P-P_water)/P;

log_coef = mu_CO2_L - log_phi_CO2 + log_gama_CO2;
m_CO2 = y_CO2*P/exp(log_coef);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
classdef spinprefsphere < spinpreference
%SPINPRESPHERE   Class for managing preferences when solving a PDE on the
%sphere with SPINSPHERE.
%
% Available preferences ([] = defaults):
%
%   Clim                      * Limits of the colorbar when 'plot' is 'movie'.
%     []                        Default is empty, i.e., automatically chosen by
%                               the code.
%
%   colormap                  * Color look-up table. See HELP/COLORMAP.
%     ['parula']
%
%   dataplot                  * Plotting options when the solution is complex.
%     ['real']
%      'imag'
%      'abs'
%
%   dealias                   * If 'on', uses the 2/3-rule to zero out high
%     ['off']                   wavenumbers.
%      'on'
%
%   grid                      * If 'on', plots longitude and latitude circles.
%     ['off']
%      'on'
%
%   iterplot                  * Plots the solution every ITERPLOT iterations of
%     [1]                       the time-stepping loop if 'plot' is 'movie'.
%
%
%   Nplot                     * Number of grid points in each direction for
%     [256]                     the movie. If Nplot>N, the data are interpolated
%                               to a finer grid.
%
%   plot                      * Plot options: 'movie' to plot a movie of the
%     ['movie']                 the solution, 'off' otherwise.
%      'off'
%
%   scheme                    * Time-stepping scheme. HELP/IMEX for the list
%     []                        of available schemes. Default is empty, i.e.,
%                               automatically chosen by the code depending on
%                               the PDE.
%
%   view                      * Viewpoint specification when 'plot' is 'movie'.
%     [-37.5 30]
%
% Construction:
%
%   PREF = SPINPREFSPHERE() creates a SPINPREF2 object with the default values.
%
%   PREF = SPINPREFSPHERE(PROP1, VALUE1, PROP2, VALUE2, ...) creates a
%   SPINPREFSPHERE object with the properties PROP1 and PROP2 set to VALUE1 and
%   VALUE2.
%
% See also SPINSPHERE.

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

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS PROPERTIES:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties ( Access = public )
        Clim                  % Limits of the colorbar (1x2*NVARS DOUBLE)
        colormap              % Color look-up table (STRING)
        grid                  % For longitude and latitude circles (STRING)
        view = [-37.5 30];    % Viewpoint of the plot (1x2 DOUBLE)
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function pref = spinprefsphere(varargin)
            if ( nargin == 0 )
                pref.colormap = 'parula';
                pref.grid = 'off';
                pref.dataplot = 'real';
                pref.dealias = 'off';
                pref.iterplot = 1;
                pref.Nplot = 128;
                pref.plot = 'movie';
            else
                pref = spinprefsphere();
                for k = 1:nargin/2
                    pref.(varargin{2*(k-1)+1}) = varargin{2*k};
                end
            end
        end
    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [newsubs, newsz] = renumber(subs, sz, range)
%RENUMBER indices for sptensor subsref
%
%  [NEWSUBS,NEWSZ] = RENUMBER(SUBS,SZ,RANGE) takes a set of
%  original subscripts SUBS with entries from a tensor of size
%  SZ. All the entries in SUBS are assumed to be within the
%  specified RANGE. These subscripts are then renumbered so that,
%  in dimension i, the numbers range from 1:numel(RANGE(i)).
%
%  See also SPTENSOR/SUBSREF
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

newsz = sz;
newsubs = subs;
for i = 1 : size(sz,2)
    if ~(ischar(range{i}) && range{i} == ':')
	if (isempty(subs))
	    newsz(i) = numel(range{i});
	else
	    [newsubs(:,i), newsz(i)] = ...
		renumberdim(subs(:,i), sz(i), range{i});
	end
    end
end

%------------------------------------------------------
function [newidx, newsz] = renumberdim(idx, sz, range)
%RENUMBERDIM helper function for RENUMBER
%  See also SPTENSOR/PRIVATE/RENUMBER

% Determine the size of the new range
newsz = numel(range);

% Create a map from the old range to the new range
map = zeros(1, sz);
for i = 1 : newsz
    map(range(i)) = i;
end

% Do the mapping
newidx = map(idx);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function gf=gradloglikGaP(x, varargin)
% f=loglikGaP(x, varargin)
% complete log likellihood of the GaP model
% Vxt = varargin{1};      %data
% sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
% alpha = varargin{3}; %parameters of the Gamma prior on the blinking
% beta = varargin{4}; %parameters of the Gamma prior on the blinking
% peval = varargin{5}; %parameters
% x(1:end-2*peval.ncomp) is Hkt

Vxt = varargin{1};      %data
sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
alpha = varargin{3}; %parameters of the Gamma prior on the blinking
beta = varargin{4}; %parameters of the Gamma prior on the blinking
peval = varargin{5}; %parameters

% % % Hkt_linear=x(1:end-peval.ncomp*2); % intensities
% % % cx=x(end-peval.ncomp*2+1:end-peval.ncomp); %x-coordinates of the centers
% % % cy=x(end-peval.ncomp+1:end); % y-coordinates of the centers
% % %
Hkt_linear=x; % intensities
cx=varargin{6};
cy=varargin{7};
% % %

sigpsf_vec=repmat(sigpsf,peval.ncomp,1); %all psfs same sigma
cxy_vec=[cx'+.5, cy'+.5];
a_vec=1./(sigpsf_vec*2*pi); % all normalised to 1

Hkt=reshape(Hkt_linear, peval.ncomp, peval.nt);
% generate PSFs from given parameters:
Wxkpix=gauss2dmultislice([peval.nx, peval.ny, peval.ncomp], cxy_vec, sigpsf_vec, a_vec);
Wxkpix=normalizePSF(Wxkpix); %normalize PSFs to 1
Wxk=reshape(Wxkpix,peval.nx*peval.ny, peval.ncomp);
% add the background component (just for muttiplication - in the gradient it should not matter):
[Wxkbg,Hktbg]=addbg(Wxk, Hkt, peval.bg);
P=Wxkbg*Hktbg; %current approximation

xxvc = lineargrad([peval.nx, peval.ny, peval.ncomp], cx, 'xx'); %linera grasdient shifted by cx
yyvc = lineargrad([peval.nx, peval.ny, peval.ncomp], cy, 'yy');

% dW/dcx:
Wxtcx=1/sigpsf^2*xxvc.*Wxk;
% dW/dcy:
Wxtcy=1/sigpsf^2*yyvc.*Wxk;

% d(log(L))/dHkt:
% % % gfHkt=(alpha-1)*1./Hkt - 1/beta + Wxk'*(Vxt./P-eye(peval.nx*peval.ny,peval.nt));
% % % gfHkt= Wxk'*(Vxt./P-eye(peval.nx*peval.ny,peval.nt));
gfHkt= Wxk'*(Vxt./P)-1; %without background
% d(log(L))/dcx:
gfcx=diag(Wxtcx'*(Vxt./P)*Hkt');
% d(log(L))/dcy:
gfcy=diag(Wxtcy'*(Vxt./P)*Hkt');

% % % gf = [reshape(gfHkt',1,peval.nt*peval.ncomp), gfcx', gfcy'];
% gf = [reshape(gfHkt(1:peval.ncomp, :),1,peval.nt*peval.ncomp)];
gf = [reshape(gfHkt,1,peval.nt*peval.ncomp)];

end
function Wnorm=normalizePSF(W)
sw=size(W);
Wr=reshape(W, sw(1)*sw(2),sw(3));
q=squeeze(sum(Wr,1));
Wrnorm=Wr./repmat(q,sw(1)*sw(2),1);
Wnorm=reshape(Wrnorm,sw(1), sw(2), sw(3));
end
function xxvc = lineargrad(sizevec, cx, dir)
switch dir
    case 'xx'
        xxp=double(xx(sizevec, 'true')); %linear function - pixels
    case 'yy'
        xxp=double(yy(sizevec, 'true')); %linear function - pixels
    otherwise
        error('Wrong dir')
end

xxv=reshape(xxp,sizevec(1)*sizevec(2),sizevec(3)); %linear function - vector
xxvc=xxv-repmat(cx,sizevec(1)*sizevec(2),1);
% xxp=double(xx([peval.nx, peval.ny, peval.ncomp], 'true')); %linear function - pixels
%yyp=double(yy([peval.nx, peval.ny, peval.ncomp], 'true'));
% xxv=reshape(xxp,peval.nx*peval.ny,peval.ncomp); %linear function - vector
%yyv=reshape(yyp,peval.nx*peval.ny,peval.ncomp);
% xxvc=xxv-repmat(cx,peval.ncomp,1);
%yyvc=yyv-repmat(cy,peval.ncomp,1);
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [ Xt,Yt,Zt ] = multiplyEuMat( EuMat, X,Y,Z )
%MULTIPLYEUMAT takes the X, Y, Z coordinates of an object and returns the
%coordinates Xt, Yt, Zt of same object rotated using a rotation matrix

Xt=X;
Yt=Y;
Zt=Z;

resvec=[1;1;1];
for i=1:numel(X)
   temp=[X(i);Y(i);Z(i)];
   resvec=EuMat*temp;
   Xt(i)=resvec(1);
   Yt(i)=resvec(2);
   Zt(i)=resvec(3);
end

author: Sable

repo: mcbench-benchmarks


text:
function I2 = im_fusion(I,wl,I1,options)
%IM_FUSION Summary of this function goes here
%   Detailed explanation goes here
% Input:
%   I - rows by cols by B hyperspectral image, in the range [0 1]
%   wl - 1 by B wavelength
%   I1 - rows1 by cols1 by b RGB image, in the range [0 255] or [0 1]
if nargin < 4
    options = [];
end

t_start = tic;

I1 = double(I1);
scale_rgb = 0;
if max(I1(:)) > 10 % in the range 0 - 255
    I1 = I1/255;
    scale_rgb = 1;
end

has_reg_params = parse_param(options,'has_reg_params',0);
bbl = parse_param(options,'bbl',true(1,size(I,3)));
show_fig = parse_param(options,'show_fig',0);
has_PSF_SRF = parse_param(options,'has_PSF_SRF',0);
regularization = parse_param(options,'regularization','LLMM'); %LMM or LLMM
sylvester_solver_mode = parse_param(options,'sylvester_solver_mode','adaptive_accurate');

%% automatically infer s and extra and check size match
[rows,cols,B] = size(I);
[rows1,cols1,b] = size(I1);

s = parse_param(options,'s', [floor(cols1/cols),floor(rows1/rows)]);
extra = parse_param(options,'extra',max(([cols1,rows1]-s.*[cols,rows])/2, [1,1]));

pad_color = 0;
if all([cols,rows].*s + extra*2 == [cols1,rows1])
    disp('Image sizes match exactly');
elseif all([cols,rows].*s == [cols1,rows1])
    % add padding for extra margin
    I1 = pad_rgb_margin(I1, extra);
    disp('Pad color image with extra margin');
    [rows1,cols1,b] = size(I1);
    pad_color = 1;
else
    error('Image sizes do not match');
end

%% decode PSF and SRF
if has_PSF_SRF
    H1 = options.H1;
    G = options.G;
    S = options.S;
    s = options.s;
    g = options.g;
    extra = options.extra;
    if scale_rgb
        H1 = H1/255;
    end
    I(:,:,~bbl) = [];
    wl(~bbl) = [];
elseif has_reg_params
    reg = options.reg_params;
    % bad bands are removed in interpret_reg_params
    [I,wl,I1,g,H1,extra,s,G,S] = interpret_reg_params(I,wl,bbl,...
        I1,reg.U2,reg.V2,reg.rho,reg.s2,reg.T2,reg.sigma2);
else
    % bad bands are removed in estimate_PSF_SRF
    [I,wl,I1,g,H1,extra,s,G,S] = estimate_PSF_SRF(I,wl,bbl,I1,s,extra,options);
end

%% adjust parameters
[rows,cols,B] = size(I);
[rows1,cols1,b] = size(I1);

N = rows * cols;
N1 = rows1 * cols1;

h0 = H1(1,:)';
H = H1(2:end,:);
Y = reshape_hsi(I);
X = reshape(I1, [N1, b]);

beta = parse_param(options,'beta',1);
% beta = max(beta, 1e-12);
%
% beta2 = parse_param(options,'beta2',0.1);
% tau = beta2/beta;

beta = beta * (b/B);
% options.tau = tau;

gamma = parse_param(options,'gamma',0.5);%weight
gamma = max(gamma,1e-4);
gamma = 1/(N*B*(1-gamma)/(N1*b*gamma) + 1);

%% Do fusion
disp('Start image fusion');

F = S*H;
X1 = X - ones(N1,1) * h0';
Z = gamma*G'*Y + (1-gamma)*X1*F';

options.regularization = regularization;
L = calc_regularization(reshape(X1,[rows1,cols1,b]), options);

options = insert_into_options(options,'g',g, ...
    'beta',beta,'gamma',gamma,'G',G,'Y',Y,'F',F,'X1',X1,'Z',Z, ...
    'extra',extra,'s',s,'L',L);

% R = optimize_gradient_descent(I, I1, options);
% R = fusion_with_couple_Local_LMM(I,I1,options);

if strcmp(regularization, 'LLMM')
    if 1 % reconstruct full reflectance data
        se_A = gamma*(G'*G) + beta*L;
        se_B = (1-gamma)*(F*F');
        se_C = Z;
        R = sylvester_krylov(se_A, se_B, se_C, sylvester_solver_mode);
    else % reconstruct abundances
        num_endm = vd(Y',5*10^-2); % M can be automatically defined, for example, by VD
        [~, indicies] = vca(Y', num_endm);
        M = Y(indicies,:);
        [U,Lambda] = eig(M*M');
        tmp = M*(1-gamma)*F*F'*M';
        [V,S] = eig((tmp + tmp')/2);
        Lambda1 = diag(sqrt(1./diag(Lambda)));
        se_A = gamma*(G'*G) + beta*L;
        se_B = Lambda1*U'*V*S*V'*U*Lambda1;
        se_B = (se_B + se_B) / 2;
        se_C = Z*M'*U*Lambda1;
        Z = sylvester_krylov(se_A, se_B, se_C, 'adaptive');
        A = Z*Lambda1*U';
        R = A*M;
    end
elseif strcmp(regularization, 'LMM')
    kappa = parse_param(options,'kappa',0.1);
    kappa = kappa * (b/B);
    options.kappa = kappa;

    R = fusion_with_LMM_constraint(I, I1, options);
end

%% pad bad bands
if isfield(options,'bbl') && ~all(options.bbl)
    R1 = nan(size(R,1),length(options.bbl));
    R1(:,options.bbl) = R;
    R = R1;
end

I2 = reshape(R, [rows1,cols1,size(R,2)]);

if pad_color
    I2 = remove_rgb_margin(I2, extra);
end

t_elapsed = toc(t_start);
disp(['Elapsed time for image fusion is ', num2str(t_elapsed)]);
end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function fem = update11nbr(fem_struct,jj)
% update11nbr takes a node with 11 elements connected to it
%  (which defines a region) and addes 4 nodes and 8 elements
%  to that region then regroups so that no node has more than
%  six elements connected to it. The new nodes and elements are sprung.
%
% NOTE:  fem_struct.nei must be a component.
%
%  Variables
%  fem_struct -- finite element structure from opnml.
%  jj -- the node number that has 11 elements connected to it.
%  fem -- the updated finite element structure.
%        Note:  Only fem.x,fem.y,fem.z,fem.e, and fem.nei
%               get updated.  fem.bnd does not need to be updated.
%
%  Usage -- fem = update11nbr(fem_struct,jj)
%
%  Name: update11nbr.m
%  Written by: Ben Holladay
%  Date: June 22,2004
%  Updated:  Aug. 27, 2004 -- Chris Massey Fixed bug in .nei for the
%                                          3rd new node.
% Last Modifed:
%         Sept. 19, 2006 -- Chris Massey, NRL Code 7322, S.S.C. MS 39529
%                      % Add test to make sure the correct nodal neighbor
%                        connectivity existed for the requested update
%                        and added a test to create the nodal neighbor list
%                        if not present.
%         Aug. 13, 2007 -- Chris Massey, NRL Code 7322
%                     Fixed a bug in the code that updates nodal neighbor
%                     list for nodes getting extra connectivity.
%         July 14, 2008 -- Chris Massey, NRL Code 7322
%               moved test for 4 neighbors to after the test for nei.
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Checks if the fem_struct has an nei. If not one is generated.
try
    nei = fem_struct.nei;
catch
    fem_struct.nei = ele2nei(fem_struct.e,fem_struct.x,fem_struct.y);
    nei = fem_struct.nei;
end

tempnc = sum(~ismember((fem_struct.nei(jj,:)),0));
if tempnc ~= 11
    error(['There are ',num2str(tempnc),' nodal neighbors for node number ',...
           num2str(jj),'. This routine only works for 11 nodal neighbors.']);
end

%Sets the intial fem_struct variables.
x = fem_struct.x;
y = fem_struct.y;
z = fem_struct.z;
enodes = fem_struct.e;
vod = [1:11,1:11];

%Adds new nodes and elements and labels the bad node, neighboring nodes,
%and neighboring elements.
[nnodes,nc] = size(nei);
nelems = size(enodes,1);
newnode = nnodes + [1:4];
newelem = nelems + [1:8];
badnode = jj;
nbrnode = nei(jj,1:11);
[nbrelem,~] = find(enodes == badnode);
inrad = max(sqrt((x(nbrnode)-x(badnode)).^2 + (y(nbrnode)-y(badnode)).^2));

%First guess of the new coordinates and bathymetry.
Bx = sum(x([(nbrnode([1:3])),badnode]))/4;
By = sum(y([(nbrnode([1:3])),badnode]))/4;
Bz = sum(z([(nbrnode([1:3])),badnode]))/4;

Ax = sum([(x([(nbrnode([10,11,1]))]))',Bx])/4;
Ay = sum([(y([(nbrnode([10,11,1]))]))',By])/4;
Az = sum([(z([(nbrnode([10,11,1]))]))',Bz])/4;

Cx = sum([(x([(nbrnode([3:6]))]))',Bx])/5;
Cy = sum([(y([(nbrnode([3:6]))]))',By])/5;
Cz = sum([(z([(nbrnode([3:6]))]))',Bz])/5;

Ex = sum(x([(nbrnode([7:10])),badnode]))/5;
Ey = sum(y([(nbrnode([7:10])),badnode]))/5;
Ez = sum(z([(nbrnode([7:10])),badnode]))/5;

Dx = sum([(x([(nbrnode([6:7]))]))',Bx,Ax,Cx,Ex])/6;
Dy = sum([(y([(nbrnode([6:7]))]))',By,Ay,Cy,Ey])/6;
Dz = sum([(z([(nbrnode([6:7]))]))',Bz,Az,Cz,Ez])/6;

x(badnode) = Ax;
y(badnode) = Ay;
z(badnode) = Az;
x(newnode(1)) = Bx;
y(newnode(1)) = By;
z(newnode(1)) = Bz;
x(newnode(2)) = Cx;
y(newnode(2)) = Cy;
z(newnode(2)) = Cz;
x(newnode(3)) = Dx;
y(newnode(3)) = Dy;
z(newnode(3)) = Dz;
x(newnode(4)) = Ex;
y(newnode(4)) = Ey;
z(newnode(4)) = Ez;

%Updates the effected elements. J variables are used to represent which
%elements are connected each node.
i = 1;
while i <= 9
    j1 = ismember((enodes((nbrelem),:)),(nbrnode(vod(i))));
    j2 = ismember((enodes((nbrelem),:)),(nbrnode(vod(i+1))));
    spb = [1,1,2,2,2,3,4,4,4];
    j = sum([j1,j2],2);
    temp = nbrelem(find (j == 2));
    if ~isempty(temp)
        enodes(temp,:) = ([nbrnode(vod(i)),nbrnode(vod(i+1)),newnode(spb(i))]);
    end
    i = i + 1;
    clear j j1 j2 temp;
end

%Addes the 8 new elements.
enodes(newelem(1),:) = [nbrnode(1),newnode(1),badnode];
enodes(newelem(2),:) = [nbrnode(3),newnode(2),newnode(1)];
enodes(newelem(3),:) = [newnode(1),newnode(2),newnode(3)];
enodes(newelem(4),:) = [newnode(1),newnode(3),badnode];
enodes(newelem(5),:) = [badnode,newnode(3),newnode(4)];
enodes(newelem(6),:) = [nbrnode(10),badnode,newnode(4)];
enodes(newelem(7),:) = [nbrnode(6),newnode(3),newnode(2)];
enodes(newelem(8),:) = [nbrnode(7),newnode(4),newnode(3)];

M(1,:) = ([Ax,Ay,Az]);
M(2,:) = ([Bx,By,Bz]);
M(3,:) = ([Cx,Cy,Cz]);
M(4,:) = ([Dx,Dy,Dz]);
M(5,:) = ([Ex,Ey,Ez]);

%Springs the region.
i = 1;
imax = 20;
stoptol = 10e-8 * inrad;
tol = stoptol + 10;
while i <= imax && tol > stoptol
    M2(2,1) = sum([(x([(nbrnode([1:3])),badnode]))',M(3,1),M(4,1)]);
    M2(2,2) = sum([(y([(nbrnode([1:3])),badnode]))',M(3,2),M(4,2)]);
    M2(2,3) = sum([(z([(nbrnode([1:3])),badnode]))',M(3,3),M(4,3)]);

    M2(1,1) = sum([(x([(nbrnode([10,11,1]))]))',M(2,1),M(4,1),M(5,1)]);
    M2(1,2) = sum([(y([(nbrnode([10,11,1]))]))',M(2,2),M(4,2),M(5,2)]);
    M2(1,3) = sum([(z([(nbrnode([10,11,1]))]))',M(2,3),M(4,3),M(5,3)]);

    M2(3,1) = sum([(x([(nbrnode([3:6]))]))',M(2,1),M(4,1)]);
    M2(3,2) = sum([(y([(nbrnode([3:6]))]))',M(2,2),M(4,2)]);
    M2(3,3) = sum([(z([(nbrnode([3:6]))]))',M(2,3),M(4,3)]);

    M2(5,1) = sum([(x([(nbrnode([7:10]))]))',M(4,1),M(1,1)]);
    M2(5,2) = sum([(y([(nbrnode([7:10]))]))',M(4,2),M(1,2)]);
    M2(5,3) = sum([(z([(nbrnode([7:10]))]))',M(4,3),M(1,3)]);

    M2(4,1) = sum([(x([(nbrnode([6:7]))]))',M(2,1),M(1,1),M(3,1),M(5,1)]);
    M2(4,2) = sum([(y([(nbrnode([6:7]))]))',M(2,2),M(1,2),M(3,2),M(5,2)]);
    M2(4,3) = sum([(z([(nbrnode([6:7]))]))',M(2,3),M(1,3),M(3,3),M(5,3)]);
    M2 = M2./ 6;

    tol = max(sqrt((M2(:,1)-M(:,1)).^2 + (M2(:,2)-M(:,2)).^2));

    M = M2;
    x([badnode,([newnode(1:4)])]) = ([M(:,1)]);
    y([badnode,([newnode(1:4)])]) = ([M(:,2)]);
    z([badnode,([newnode(1:4)])]) = ([M(:,3)]);
    i = i + 1;
end

%Creates the output structure.
fem = fem_struct;
fem.e = enodes;
fem.x = x;
fem.y = y;
fem.z = z;
fem.nei = nei;

%Update the connectivity for the bad node and adds the connections for
%the new nodes.
fem.nei(badnode,1:nc) = ([newnode([1,3,4]),nbrnode([10,11,1]),zeros(1,nc-6)]);
fem.nei(newnode(1),1:nc) = ([badnode,nbrnode([1:3]),newnode([2,3]),zeros(1,nc-6)]);
fem.nei(newnode(2),1:nc) = ([newnode(1),nbrnode([3:6]),newnode(3),zeros(1,nc-6)]);
fem.nei(newnode(3),1:nc) = ([badnode,newnode([1,2]),nbrnode([6:7]),newnode(4),zeros(1,nc-6)]); %TCM 8/27/04
fem.nei(newnode(4),1:nc) = ([nbrnode([7:10]),badnode,newnode(3),zeros(1,nc-6)]);

%Updates the nei for neighbor nodes. These nodes get an additional
%connector.
spb = ([1,3,6,7,10]);
for i = 1:5
    addconn = ([badnode,newnode([1:4]),badnode]);
    tmp = nei(nbrnode(spb(i)),:);
	ij = find(tmp == badnode);
    ik = min((find(tmp == 0)) + 1);
    if isempty(ik)  % TCM 08/13/2007 -- Added this check in case ik was empty.
       ik = size(nei,2)+1;
    end
    fem.nei(nbrnode(spb(i)),1:ik) = ([tmp(1:(ij-1)),addconn(i+1),addconn(i),tmp((ij+1):(ik-1))]);
end

%Updates the nei for neighbor nodes. These nodes get different connector.
spn = ([2,4,5,8,9]);
for i = 1:5
    addnode = ([newnode(1),newnode(2),newnode(2),newnode(4),newnode(4)]);
    tmp = nei(nbrnode(spn(i)),:);
    ij = find(tmp == badnode);
    fem.nei(nbrnode(spn(i)),1:length(tmp)) = ([tmp(1:(ij-1)),addnode(i),tmp((ij+1):end)]);
end

%Determine the triqual for the final form to see if has very low quality
%elements.
tmp1 = (nbrelem);
tempL3=(fem.x(fem.e(tmp1,1))-fem.x(fem.e(tmp1,2))).^2+(fem.y(fem.e(tmp1,1))-...
   fem.y(fem.e(tmp1,2))).^2;
tempL1=(fem.x(fem.e(tmp1,2))-fem.x(fem.e(tmp1,3))).^2+(fem.y(fem.e(tmp1,2))-...
   fem.y(fem.e(tmp1,3))).^2;
tempL2=(fem.x(fem.e(tmp1,3))-fem.x(fem.e(tmp1,1))).^2+(fem.y(fem.e(tmp1,3))-...
   fem.y(fem.e(tmp1,1))).^2;
tempLs = ([tempL1 + tempL2 + tempL3]);
xnodes = fem.x(fem.e(tmp1,:));
ynodes = fem.y(fem.e(tmp1,:));
temparea = 0.5*(xnodes(:,1).*(ynodes(:,2)-ynodes(:,3))+xnodes(:,2).*...
   (ynodes(:,3)-ynodes(:,1))+xnodes(:,3).*(ynodes(:,1)-ynodes(:,2)));
fem.ar(tmp1) = temparea;
tempq = (4 * sqrt(3) * temparea) ./ tempLs;
clear tempL3 tempL2 tempL1 tempLs xnodes ynodes temparea;

fem1 = fem;
%Identify very low quality elements and attempts to use a line swap to fix.
poor = find(tempq < .1);
nflag = 1;good = 1;
if ~isempty(poor)
   for it = 1:length(poor)
      je = tmp1(poor(it));
      nodes = fem.e(je,:);
      nflag = 0;

      %Determine the angles to see if an edge can be fliped
      a2 = (x(enodes(je,3))-x(enodes(je,2))).^2+(y(enodes(je,3))-y(enodes(je,2))).^2;
      b2 = (x(enodes(je,1))-x(enodes(je,3))).^2+(y(enodes(je,1))-y(enodes(je,3))).^2;
      c2 = (x(enodes(je,2))-x(enodes(je,1))).^2+(y(enodes(je,2))-y(enodes(je,1))).^2;
      A = (180/pi)*acos((b2+c2-a2)./(2*sqrt(b2).*sqrt(c2)));
      B = (180/pi)*acos((c2+a2-b2)./(2*sqrt(c2).*sqrt(a2)));
      C = (180/pi)*acos((a2+b2-c2)./(2*sqrt(a2).*sqrt(b2)));
      [test,ind] = max([A,B,C]);
      if test > 160
         %Find the element numbers to flip the edge.
         nogood = nodes(ind);
         swap = setdiff(nodes,nogood);
         [temp1,tc] = find(enodes == swap(1) | enodes == swap(2));
         [b,j,k] = unique(temp1);
         temp2 = setdiff(1:length(temp1),j);
         swap = temp1(temp2);
         try
            fem1 = line_swap(fem,swap(1),swap(2));
         catch
            fem1 = fem;
         end

         %Spring the new mesh after the line swap.
         for itt = 1:2
             temp = find(fem1.nei(newnode(1),:) ~= 0);
             tempnei = fem1.nei(newnode(1),temp);
             fem1.x(newnode(1)) = mean(fem1.x(tempnei));
             fem1.y(newnode(1)) = mean(fem1.y(tempnei));
             fem1.z(newnode(1)) = mean(fem1.z(tempnei));
             temp = find(fem1.nei(newnode(2),:) ~= 0);
             tempnei = fem1.nei(newnode(2),temp);
             fem1.x(newnode(2)) = mean(fem1.x(tempnei));
             fem1.y(newnode(2)) = mean(fem1.y(tempnei));
             fem1.z(newnode(2)) = mean(fem1.z(tempnei));
             temp = find(fem1.nei(newnode(3),:) ~= 0);
             tempnei = fem1.nei(newnode(3),temp);
             fem1.x(newnode(3)) = mean(fem1.x(tempnei));
             fem1.y(newnode(3)) = mean(fem1.y(tempnei));
             fem1.z(newnode(3)) = mean(fem1.z(tempnei));
             temp = find(fem1.nei(newnode(4),:) ~= 0);
             tempnei = fem1.nei(newnode(4),temp);
             fem1.x(newnode(4)) = mean(fem1.x(tempnei));
             fem1.y(newnode(4)) = mean(fem1.y(tempnei));
             fem1.z(newnode(4)) = mean(fem1.z(tempnei));
             temp = find(fem1.nei(badnode,:) ~= 0);
             tempnei = fem1.nei(badnode,temp);
             fem1.x(badnode) = mean(fem1.x(tempnei));
             fem1.y(badnode) = mean(fem1.y(tempnei));
             fem1.z(badnode) = mean(fem1.z(tempnei));
         end

         %Use triqual to determine if the new mesh is better quality.
         tmp1 = nbrelem;
         tempL3=(fem1.x(fem1.e(tmp1,1))-fem1.x(fem1.e(tmp1,2))).^2+(fem1.y(fem1.e(tmp1,1))-...
            fem1.y(fem1.e(tmp1,2))).^2;
         tempL1=(fem1.x(fem1.e(tmp1,2))-fem1.x(fem1.e(tmp1,3))).^2+(fem1.y(fem1.e(tmp1,2))-...
            fem1.y(fem1.e(tmp1,3))).^2;
         tempL2=(fem1.x(fem1.e(tmp1,3))-fem1.x(fem1.e(tmp1,1))).^2+(fem1.y(fem1.e(tmp1,3))-...
            fem1.y(fem1.e(tmp1,1))).^2;
         tempLs = ([tempL1 + tempL2 + tempL3]);
         xnodes = fem1.x(fem1.e(tmp1,:));
         ynodes = fem1.y(fem1.e(tmp1,:));
         temparea = 0.5*(xnodes(:,1).*(ynodes(:,2)-ynodes(:,3))+xnodes(:,2).*...
            (ynodes(:,3)-ynodes(:,1))+xnodes(:,3).*(ynodes(:,1)-ynodes(:,2)));
         fem1.ar(tmp1) = temparea;
         tempq = (4 * sqrt(3) * temparea) ./ tempLs;
         clear tempL3 tempL2 tempL1 tempLs xnodes ynodes temparea;
         if min(tempq) > .4
            fem = fem1;
            nflag = 1;
         else
            nflag = 0;
            good = 6;
         end
      end
   end
end

% If problems with new mesh, then retriangulate

%Correct output if line swap failed to run
if sum(nflag) == 0
   fem = patch_update(fem_struct,fem1,jj);
   good = 6;
end

% Correct output if invalid elements were created.
bad = find(fem.ar < 0);
if ~isempty(bad)
    fem = patch_update(fem_struct,fem1,jj);
    good = 6;
end

%Display message if invalid elements where created.
if good == 6
   disp('The nodally updated mesh contained invalid or badly conditioned');
   disp('elements, therefore the patch was retriangulated which should');
   disp('reduce the connecitivity but is not guaranteed to do so.');
   disp(' ');
end

return

author: CHLNDDEV