headplot.m

来自「含有多种ICA算法的eeglab工具箱」· M 代码 · 共 797 行 · 第 1/2 页

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if strcmpi(g.orilocs, 'off')
        fprintf('Computing electrode locations on head...\n');
        for i=1:length(Xe)
            elect = [Xe(i) Ye(i) Ze(i)];
            dists = distance(elect,spherePOS');
            [S,I] = sort(dists);
            npoints = I(1:3);
            diffe = newPOS(npoints,:)-spherePOS(npoints,:);
            newElect(i,:) = elect+mean(diffe)*ElectDFac;
            %if Ze(i) < 0               % Plot superior electrodes only.
            %    newElect(i,:) = [0 0 0]; % Mark lower electrodes  as having
            %end                        % an electrode position not to be plotted
        end
    else 
        fprintf('Using original electrode locations on head...\n');
        newElect(:,1) = -cell2mat( tmpY(indices) )';
        newElect(:,2) =  cell2mat( tmpX(indices) )';
        newElect(:,3) =  cell2mat( tmpZ(indices) )';        
    end;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Calculate g(x) for sphere mesh vertices
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fprintf('Computing %d vertices. Should take a while (see wait bar)\n',...
                      length(x))
    fprintf('            but doesnt have to be done again for this montage...\n');
    hwb = waitbar(0,'Computing spline file (percent done)...');

    %figure; plot3(x (:),  y(:),  z(:), '.' );
    %figure; plot3(Xe(:), Ye(:), Ze(:), '.' );

    hwbend = length(x);
    for j = 1:length(x)
      % fprintf('%d ',j)
      X = x(j);
      Y = y(j);
      Z = z(j);
      ei = onemat-sqrt((X*onemat-Xe).^2 + (Y*onemat-Ye).^2 + (Z*onemat-Ze).^2); 
                                  % default was /2, no sqrt
      for i = 1:length(ei)
        gx(j,i) = calcgx(ei(i));  
      end
      waitbar(j/hwbend)
    end
    close(hwb)

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Save spline file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    comment = g.comment;
    save(spline_file, '-mat', 'Xe', 'Ye', 'Ze', 'G', 'gx', 'newElect', 'ElectrodeNames', 'comment');
    tmpinfo = dir(spline_file);
    fprintf('Saving (%dk) file %s\n',round(tmpinfo.bytes/1000), spline_file);
    return

  elseif strcmp(values,'example') | strcmp(values,'demo')
%
%%%%%%%%%%%%%%%%%% Show an example electrode angles file  %%%%%%%%%%%%%%%%%%%%%%%%
%
       fprintf(['\nExample of a headplot() electrode angles file (spherical coords.)\n',...
               'Fields:  chan_num cor_deg horiz_deg channel_name\n\n',...
               '           1        -90     -72        Fp1.\n',...
               '           2         90      72        Fp2.\n',...
               '           3        -62     -57        F3..\n',...
               '           4         62      57        F4..\n',...
               '           5        -45       0        C3..\n',...
               '           6         45       0        C4..\n',...
               '           7       -118       2        A1..\n',...
               '           8        118      -2        A2..\n',...
               '           9        -62      57        P3..\n',...
               '           10        62     -57        P4..\n',...
               '           11       -90      72        O1..\n',...
               '           12        90     -72        O2..\n',...
               '           13       -90     -36        F7..\n',...
               '           14        90      36        F8..\n',...
               '           15       -90       0        T3..\n',...
               '           16        90       0        T4..\n',...
               '           17       -90      36        T5..\n',...
               '           18        90     -36        T6..\n',...
               '           19        45      90        Fz..\n',...
               '           20         0       0        Cz..\n',...
               '           21        45     -90        Pz..\n',...
             '\nA 90 deg coronal rotation points to right ear, -90 to left.\n' ,...
               'A positive horizontal rotation is counterclockwise from above.\n',...
               'Use pol2sph() to convert from topoplot() format to spherical.\n',...
               'Channel names should have 4 chars (. = space).\n\n\n']);
       return
  elseif strcmp(values,'cartesian') 
%
%%%%%%%%%%%%%%%%%% Show an example cartesian electrode file  %%%%%%%%%%%%%%%%%%%
%
       fprintf(['\nExample of a headplot() electrode location file (cartesian coords.)\n',...
               'Fields:  chan_num  x        y        z        channel_name\n\n',...
               '           1       0.4528   0.8888  -0.0694          Fp1.\n',...
               'Channel names should have 4 chars (. = space).\n\n\n']);
       return
  else
    fprintf('headplot(): Unknown first argument (%s).\n',values)
    help headplot
  end
else
%
%%%%%%%%%%%%%%%%%%%%%%%%%% Make the plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   
   if nargin < 2
       help headplot
       return
   end
   spline_file = arg1;
   
   g = finputcheck( varargin, { ...
       'cbar'       'real'   [0 Inf]         []; % Colorbar value must be 0 or axis handle.'
       'lighting'   'string' { 'on' 'off' }  'on';
       'verbose'    'string' { 'on' 'off' }  'on';
       'maplimits'  { 'string' 'real' }  []  'absmax'; 
       'title'      'string' []              '';
       'lights'     'real'   []              DEFAULT_LIGHTS;
       'view'       'real'   []              [143 18];
       'colormap'   'real'   []              jet(64);
       'meshfile'   'string' []              DEFAULT_MESH;
       'electrodes' 'string' { 'on' 'off' }  'on';            
       'orilocs'    { 'string' 'struct' } [] '';            
       'labels'     'integer' [0 1 2]        0 }, 'headplot');
   if isstr(g) error(g); end;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Open head mesh and electrode spline files
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if ~exist(spline_file)
       error(sprintf('headplot(): spline_file "%s" not found. Run headplot in "setup" mode\n',...
           spline_file));
  end
  load(spline_file, '-mat');
  enum = length(values);
  if enum ~= length(Xe)
	  close;
	  error('headplot(): Number of values in spline file should equal number of electrodes')
  end
  % change electrode if necessary
  % -----------------------------
  if ~isempty(g.orilocs)
      eloc_file = readlocs( g.orilocs );
      fprintf('Using original electrode locations on head...\n');
      indices = find(~cellfun('isempty', { eloc_file.X } ));
      newElect(:,1) = -cell2mat( { eloc_file(indices).Y } )';
      newElect(:,2) =  cell2mat( { eloc_file(indices).X } )';
      newElect(:,3) =  cell2mat( { eloc_file(indices).Z } )';        
  end;
  
  % load mesh file
  % --------------
  if ~exist(g.meshfile)
      error(sprintf('headplot(): mesh file "%s" not found\n',g.meshfile));
  end
  try, load(g.meshfile,'-mat');
  catch,
      POS  = load('mheadnewpos.txt', '-ascii');
      TRI1 = load('mheadnewtri1.txt', '-ascii'); % upper head
      % TRI2 = load('mheadnewtri2.txt', '-ascii'); % lower head
      % index1 = load('mheadnewindex1.txt', '-ascii');
      center = load('mheadnewcenter.txt', '-ascii');
  end;
  if exist('index1') ~= 1, index1 = sort(unique(TRI1(:))); end;
  if exist('TRI2')   ~= 1, TRI2 = []; end;
  if exist('NORM')   ~= 1, NORM = []; end;
  if exist('TRI1')   ~= 1, error('Variable ''TRI1'' not defined in mesh file'); end;
  if exist('POS')    ~= 1, error('Variable ''POS'' not defined in mesh file'); end;
  if exist('center') ~= 1, center = [0 0 0]; disp('Using [0 0 0] for center of head mesh'); end;

  %%%%%%%%%%%%%%%%%%%%%%%%%%
  % Perform interpolation
  %%%%%%%%%%%%%%%%%%%%%%%%%%

  meanval = mean(values); values = values - meanval; % make mean zero
  onemat = ones(enum,1);
  lamd = 0.1;
  C = pinv([(G + lamd);ones(1,enum)]) * [values(:);0]; % fixing division error
  P = zeros(1,size(gx,1));
  for j = 1:size(gx,1)
    P(j) = dot(C,gx(j,:));
  end
  P = P + meanval;
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%
  % Plot surfaces
  %%%%%%%%%%%%%%%%%%%%%%%%%%
  cla % clear axis
  HeadAxes = gca;

  W = zeros(1,size(POS,1));
  m = size(g.colormap,1);
  if size(g.maplimits) == [1,2]
      amin = g.maplimits(1);
      amax = g.maplimits(2);
  elseif strcmp(g.maplimits,'maxmin') | strcmp(g.maplimits,'minmax')
      amin = min(min(abs(P)))*1.02; % 2% shrinkage keeps within color bounds
      amax = max(max(abs(P)))*1.02; 
  elseif strcmp(g.maplimits,'absmax')
      amin = min(min(abs(P)))*1.02; % 2% shrinkage keeps within color bounds
      amax = max(max(abs(P)))*1.02; 
      amax = max(-amin, amax);
      amin = -amax;
      %amin = -max(max(abs(P)))*1.02; % 2% shrinkage keeps within color bounds
      %amax = -amin;    
  end 
  
  idx = min(m,round((m-1)*(P-amin)/(amax-amin))+1); % get colormap indices
  %subplot(1,2,1); hist(P(:));
  %idx = round((m-1)*P/(amax-amin))+m/2;
  %idx = max(1,min(m,idx)); % get colormap indices
  %subplot(1,2,2); hist(idx(:)); 
  %return;

  W(index1) = idx;
  colormap(g.colormap)
  p1 = patch('Vertices',POS,'Faces',TRI1,'FaceVertexCdata',W(:),...
      'FaceColor','interp', 'cdatamapping', 'direct', 'tag', 'mesh');    %%%%%%%%% Plot scalp map %%%%%%%%%
  if exist('NORM') == 1 & ~isempty(NORM)
      set(p1, 'vertexnormals', NORM);
  end;
  
  if ~isempty(TRI2)
      FCmap = [g.colormap; g.colormap(end,:); FaceColor; FaceColor; FaceColor];
      colormap(FCmap)
      W = ones(1,size(POS,1))*(m+2);
      p2 = patch('Vertices',POS,'Faces',TRI2,'FaceColor','interp',...
                 'FaceVertexCdata',W(:)); %%%%%%%% Plot face and lower head %%%%%%
  else 
      p2 = [];
  end;

  axis([-125 125 -125 125 -125 125])
  axis off % hide axis frame
  
  %%%%%%%%%%%%%%%%%%%%%%%%%
  % Draw colorbar - Note: uses enhanced cbar() function by Colin Humphries
  %%%%%%%%%%%%%%%%%%%%%%%%%      
  if ~isempty(g.cbar)
      BACKCOLOR = get(gcf,'Color');
      if g.cbar == 0       
          ColorbarHandle = cbar(0,3,[amin amax]); 
      else
          ColorbarHandle = cbar(g.cbar,3,[amin amax]); 
      end
      pos = get(ColorbarHandle,'position');  % move left & shrink to match head size
      set(ColorbarHandle,'position',[pos(1)-.05 pos(2)+0.13 pos(3)*0.7 pos(4)-0.26]);
  end
  axes(HeadAxes);

  %%%%%%%%%%%%%%%%%%%%%%%%%
  % Turn on lights
  %%%%%%%%%%%%%%%%%%%%%%%%%

  if strcmp(g.lighting,'on')
    set([p1 p2],'EdgeColor','none')
    
    for i = 1:size(g.lights,1)
      hl(i) = light('Position',g.lights(i,:),'Color',[1 1 1],...
      'Style','infinite');
    end
    if ~isempty(p2)
        set(p2,'DiffuseStrength',.6,'SpecularStrength',0,...
               'AmbientStrength',.4,'SpecularExponent',5)
    end;
    set(p1,'DiffuseStrength',.6,'SpecularStrength',0,...
    'AmbientStrength',.3,'SpecularExponent',5)
    lighting phong  % all this gives a matte reflectance
  end  

  %%%%%%%%%%%%%%%%%%%%%%%%%
  % Set viewpoint
  %%%%%%%%%%%%%%%%%%%%%%%%%

  if isstr(g.view)
    switch lower(g.view)
      case {'front','f'}
        view(-180,30)
      case {'back','b'}
        view(0,30)
      case {'left','l'}
        view(-90,30)
      case {'right','r'}
        view(90,30)
      case {'frontright','fr'}
        view(135,30)
      case {'backright','br'}
        view(45,30)
      case {'frontleft','fl'}
        view(-135,30)
      case {'backleft','bl'}
        view(-45,30)
      case 'top'
        view(0,90)
      case 'bottom'    % undocumented option!
        view(0,-90)
        Lights = [-125 125 80;   ...
                   125 125 80;   ...
                   125 -125 125; ...
                  -125 -125 125; ...
                   0 10 -80]; % add light from below!
      otherwise
        close; error(['headplot(): Invalid View value %s',g.view])
    end
  else
      if ~isstr(g.view)
          [h,a] = size(g.view);
          if h~= 1 | a~=2
              close; error('headplot(): View matrix size must be (1,2).')
          end
      end
      view(g.view)   % set camera viewpoint
  end
  
  if strcmp(g.electrodes,'on') % plot the electrode locations
   if exist('newElect')
    newNames = newElect*NamesDFac; % Calculate electrode label positions
    for i = 1:size(newElect,1)
        if newElect(i,:) ~= [0 0 0]  % plot radial lines to electrode sites
            line([newElect(i,1) HeadCenter(1)],[newElect(i,2) HeadCenter(2)],...
                 [newElect(i,3) HeadCenter(3)],'color',MarkerColor,'linewidth',1);
            
            if g.labels == 1        % plot electrode numbers
                t=text(newNames(i,1),newNames(i,2),newNames(i,3),int2str(i)); 
                set(t,'Color',NamesColor,'FontSize',NamesSize,'FontWeight','bold',...
                      'HorizontalAlignment','center');
                
            elseif g.labels == 2   % plot electrode names
                if exist('ElectrodeNames')
                    name = sprintf('%s',ElectrodeNames(i,:));
                    t=text(newNames(i,1),newNames(i,2),newNames(i,3),name);
                    set(t,'Color',NamesColor,'FontSize',NamesSize,'FontWeight','bold',...
                          'HorizontalAlignment','center'); 
                else
                    fprintf('Variable ElectrodeNames not read from spline file.\n');
                end
            else               % plot electrode markers
                line(newElect(:,1),newElect(:,2),newElect(:,3),'marker',...
                     '.','markersize',20,'color',MarkerColor,'linestyle','none');
            end
        end
    end
   else
       fprintf('Variable newElect not read from spline file.\n');
   end
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Turn on rotate3d, allowing rotation of the plot using the mouse
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if strcmp(g.verbose,'on')
    rotate3d on;   % Allow 3-D rotation of the plot by dragging the
  else             % left mouse button while cursor is on the plot
    rotate3d off
  end              
  % Make axis square
  if sqaxis
    axis image    % keep the head proportions human and as large as possible
  end
  % Add a plot title
  if ~isempty(g.title);
    % title(['\n' g.title],'fontsize',title_font);
    title([g.title],'fontsize',title_font); % Note: \n not interpreted by matlab-5.2
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  calcgx() - function used in 'setup'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [out] = calcgx(in)

out = 0;
m = 4;       % 4th degree Legendre polynomial
for n = 1:7  % compute 7 terms
  L = legendre(n,in);
    out = out + ((2*n+1)/(n^m*(n+1)^m))*L(1);
end
out = out/(4*pi);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  distance() - function used in 'setup'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [out] = distance(w,p)
% w is a matrix of row vectors
% p is a matrix of column vectors

l1 = size(w,1);
l2 = size(p,2);
out = zeros(l1,l2);

for i = 1:l1
  x = w(i,:)'*ones(1,l2);
  out(i,:) = sum((x-p).^2).^.5;
end