% HEADPLOT - plot a spherically-splined EEG field map on a semi-realistic
% 3-D head model. Can 3-D rotate the head image using the left
% mouse button.
% Example:
% >> headplot example % show an example spherical 'eloc_angles' file
% >> headplot cartesian % show an example cartesian 'eloc_angles' file
%
% Setup usage (do only once for each scalp montage):
%
% >> headplot('setup', elocs, splinefile, 'Param','Value',...);
% %
% % NOTE: previous call format below is still supported
% % >> headplot('setup', elocs, splinefile, comment, type);
%
% Required Setup-mode Inputs:
%
% elocs - file of electrode locations (compatible with READLOCS),
% or EEG.chanlocs channel location structure. If the channel
% file extension is not standard, use READLOCS to load the
% data file, e.g.
% >> headplot('setup', ...
% readlocs('myfile.xxx', 'filetype', 'besa'),...
% 'splinefile');
% splinefile - name of spline file to save spline info into. It is saved as a
% *.mat file and should be given the extension .spl .
%
% Optional Setup-mode Inputs:
%
% 'meshfile' - ['string' or structure] Matlab files containing a mesh. The
% mesh may be of different formats. It may be a Dipfit mesh as
% defined in the file standard_vol.mat. It may contain
% a structure with the fields 'vertices' and 'faces' or it
% it may contain a structure with at least two fields:
% POS - 3-D positions of vertices:
% x=left-right; y=back-front; z=up-down
% TRI1 - faces on which the scalp map should be computed
% plus possible optional fields are:
% center (optional) - 3-D center of head mesh
% TRI2 (optional) - faces in skin color
% NORM (optional) - normal for each vertex (better shading)
% 'orilocs' - ['off'|'on'] use original electrode locations on the head
% {default: 'off'} (extrapolated to spherical). Note that these
% electrode locations must be coregisted with the head mesh.
% 'transform' - [real array] Talairach-model transformation matrix to co-register
% the electrode locations with the head mesh:
% [shiftX shiftY shiftZ pitch roll yaw scaleX scaleY scaleZ]
% The transform is applied in the order shift(rotate(scale(elocs)))
% by the dipfit2.* plugin function traditionaldipfit.m
% This array is returned by COREGISTER.
% 'plotmeshonly' - [string] plot only mesh and electrode positions. Options are
% 'head' to plot the standard head mesh; 'sphere' to plot the
% texture of the head on a sphere; 'off' not to plot anything.
% {default: 'off'}
% 'comment' - ['string'] optional string containing comments for spline file
% {default: []}
%
% Standard-mode Usage thereafter:
%
% >> headplot(values,'spline_file','Param','Value',...)
%
% Required Standard-mode Inputs:
%
% values - vector containing a data value at each electrode position
% 'spline_file' - spline filename, computed and saved in 'setup' mode (above)
%
% Optional Standard-mode Inputs:
%
% 'meshfile' - [string] mesh file name. See file content in the setup-mode
% description above. {default: the EEGLAB head template file}.
% 'electrodes' - ['on'|'off'] -> show electrode positions {default 'on'}
% 'title' - Plot title {default: none}
% 'labels' - 2 -> plot stored electrode labels;
% 1 -> plot channel numbers; 0 -> no labels {default 0}
% 'cbar' - 0 -> Plot colorbar {default: no colorbar}
% Note: standard jet colormap) red = +;blue = -;green=0
% h -> Colorbar axis handle (to specify headplot location)
% 'view' - Camera viewpoint in deg. [azimuth elevation]
% 'back'|'b'=[ 0 30]; 'front'|'f'=[180 30]
% 'left'|'l'=[-90 30]; 'right'|'r'=[ 90 30];
% 'frontleft'|'bl','backright'|'br', etc.,
% 'top'=[0 90], Can rotate with mouse {default [143 18]}
% 'maplimits' - 'absmax' -> make limits +/- the absolute-max
% 'maxmin' -> scale to data range
% [min,max] -> user-definined values
% {default = 'absmax'}
% 'lights' - (3,N) matrix whose rows give [x y z] pos. of each of
% N lights {default: four lights at corners}
% 'electrode3d' - ['on'|'off'] plot electrodes in 3-D. Default is 'off'.
% 'lighting' - 'off' = show wire frame head {default 'on'}
% 'material' - [see material function] {default 'dull'}
% 'colormap' - 3-column colormap matrix {default: turbo(64)}
% 'verbose' - 'off' -> no msgs, no rotate3d {default: 'on'}
% 'orilocs' - [channel structure or channel file name] Use original
% channel locations instead of the one extrapolated from
% spherical locations. Note that if you use 'orilocs'
% during setup, this is not necessary here since the
% original channel location have already been saved.
% This option might be useful to show more channels than
% the ones actually used for interpolating (e.g., fiducials).
% 'transform' - [real array] homogeneous transformation matrix to apply
% to the original locations ('orilocs') before plotting them.
%
% Note: if an error is generated, HEADPLOT may close the current figure
%
% Authors: Arnaud Delorme, Colin Humphries, Scott Makeig, SCCN/INC/UCSD,
% La Jolla, 1998-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) Arnaud Delorme, Colin Humphries and Scott Makeig,
% CNL / Salk Institute, Feb. 1998
%
% Spherical spline method: Perrin et al. (1989) Electroenceph clin Neurophys
%
% 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.
% 12-12-98 changed electrode label lines to MarkerColor -sm
% 12-12-98 added colorbar option -sm (still graphically marred by tan rect.)
% 12-13-98 implemented colorbar option using enhanced cbar -sm
% 12-13-98 implemented 'setup' comment option -sm
% 03-20-00 added cartesian electrode locations option -sm
% 07-14-00 fixed line in CALGX -sm from -ch
% 03-23-01 documented 'cartesian' locfile option -sm
% 01-25-02 reformated help & license, added links -ad
% 03-21-02 added readlocs and the use of eloc input structure -ad
function [HeadAxes, ColorbarHandle] = headplot(values, arg1, varargin)
if nargin < 1
help headplot
return
end
%%%%%%%%%%%%%%%%%%%%%%%%%% Set Defaults %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
icadefs % load definitions
set(gca,'Color',BACKCOLOR);
DEFAULT_MESH = ['mheadnew.mat']; % upper head model file (987K)
DEFAULT_TRANSFORM = [0 -5 0 -0.1000 0 -1.5700 1040 800 950]; % stretching in different dimensions
DEFAULT_TRANSFORM = [0 -10 0 -0.1000 0 -1.600 1100 1100 1100]; % keep spherical shape.
%DEFAULT_MESH = '/home/arno/matlab/juliehiresmesh.mat';
%DEFAULT_MESH = ['/home/scott/matlab/old' '/newupper.mat']; % whole head model file (183K)
DEFAULT_LIGHTS = [-125 125 80; ...
125 125 80; ...
125 -125 125; ...
-125 -125 125]; % default lights at four corners
HeadCenter = [0 0 30];
FaceColor = [.8 .55 .35]*1.1; % ~= ruddy Caucasian - pick your complexion!
MAX_ELECTRODES = 1024;
ElectDFac = 1.06; % plot electrode marker dots out from head surface
plotelecopt.NamesDFac = 1.05; % plot electrode names/numbers out from markers
plotelecopt.NamesColor = 'k'; % 'r';
plotelecopt.NamesSize = 10; % FontSize for electrode names
plotelecopt.MarkerColor= [0 0 0];
plotelecopt.electrodes3d = 'off';
sqaxis = 1; % if non-zero, make head proportions anatomical
title_font = 18;
if ischar(values)
values = lower(values);
if strcmp(values,'setup')
%
%%%%%%%%%%%%%%%%%%% Perform splining file setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if nargin < 3
help headplot;
return;
end
eloc_file = arg1;
spline_file = varargin{1};
g = finputcheck(varargin(2:end), { 'orilocs' 'string' { 'on','off' } 'off';
'plotmeshonly' 'string' { 'head','off','sphere' } 'off';
'meshfile' { 'string','struct' } [] DEFAULT_MESH;
'chaninfo' 'struct' [] struct([]);
'plotchans' 'integer' [] [];
'ica' 'string' { 'on','off' } 'off';
'transform' 'real' [] DEFAULT_TRANSFORM;
'comment' 'string' [] '' }, 'headplot', 'ignore');
if ischar(g)
error(g);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Open electrode file
%%%%%%%%%%%%%%%%%%%%%%%%%%%
[eloc_file, labels, Th, Rd, indices] = readlocs(eloc_file);
indices = find(~cellfun('isempty', { eloc_file.X }));
% channels to plot
% ----------------
if isempty(g.plotchans), g.plotchans = [1:length(eloc_file)]; end
if ~isfield(g.chaninfo, 'nosedir'), g.chaninfo(1).nosedir = '+x'; end
indices = intersect_bc(g.plotchans, indices);
% if ICA select subset of channels if necessary
% ---------------------------------------------
if ~isfield(g.chaninfo, 'icachansind'), g.chaninfo(1).icachansind = 1:length(eloc_file); end
if strcmpi(g.ica, 'on'),
rmchans2 = setdiff_bc( g.chaninfo.icachansind, indices ); % channels to remove (non-plotted)
newinds = 1:length(g.chaninfo.icachansind);
allrm = [];
% remove non-plotted channels from indices
for index = 1:length(rmchans2)
chanind = find(g.chaninfo.icachansind == rmchans2(index));
allrm = [ allrm chanind ];
end
newinds(allrm) = [];
indices = newinds;
eloc_file = eloc_file(g.chaninfo.icachansind);
end
fprintf('Headplot: using existing XYZ coordinates\n');
ElectrodeNames = strvcat({ eloc_file.labels });
ElectrodeNames = ElectrodeNames(indices,:);
Xeori = [ eloc_file(indices).X ]';
Yeori = [ eloc_file(indices).Y ]';
Zeori = [ eloc_file(indices).Z ]';
[newPOS POS TRI1 TRI2 NORM index1 center] = getMeshData(g.meshfile);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rotate channel coordinates if necessary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmpi(lower(g.chaninfo.nosedir), '+x')
rotate = 0;
else
if strcmpi(lower(g.chaninfo.nosedir), '+y')
rotate = 3*pi/2;
elseif strcmpi(lower(g.chaninfo.nosedir), '-x')
rotate = pi;
else rotate = pi/2;
end
allcoords = (Yeori + Xeori*sqrt(-1))*exp(sqrt(-1)*rotate);
Xeori = imag(allcoords);
Yeori = real(allcoords);
end
newcoords = [ Xeori Yeori Zeori ];
%newcoords = transformcoords( [ Xe Ye Ze ], [0 -pi/16 -1.57], 100, -[6 0 46]);
%newcoords = transformcoords( [ Xeori Yeori Zeori ], g.transform(4:6), g.transform(7:9), g.transform(1:3));
% same performed below with homogeneous transformation matrix
transmat = traditionaldipfit( g.transform ); % arno
newcoords = transmat*[ newcoords ones(size(newcoords,1),1)]';
newcoords = newcoords(1:3,:)';
% original center was [6 0 16] but the center of the sphere is [0 0 30]
% which compensate (see variable Headcenter)
%newcoords = transformcoords( [ Xe Ye Ze ], -[0 0 -pi/6]);
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% normalize with respect to head center
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
newcoordsnorm = newcoords - ones(size(newcoords,1),1)*HeadCenter;
tmpnorm = sqrt(sum(newcoordsnorm.^2,2));
Xe = newcoordsnorm(:,1)./tmpnorm;
Ye = newcoordsnorm(:,2)./tmpnorm;
Ze = newcoordsnorm(:,3)./tmpnorm;
%plotchans3d([ Xe Ye Ze], cellstr(ElectrodeNames)); return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate g(x) for electrodes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmpi(g.plotmeshonly, 'off')
fprintf('Setting up splining matrix.\n');
enum = length(Xe);
onemat = ones(enum,1);
G = zeros(enum,enum);
for i = 1:enum
ei = onemat-sqrt((Xe(i)*onemat-Xe).^2 + (Ye(i)*onemat-Ye).^2 + ...
(Ze(i)*onemat-Ze).^2); % default was /2 and no sqrt
gx = zeros(1,enum);
for j = 1:enum
gx(j) = calcgx(ei(j));
end
G(i,:) = gx;
end
end
fprintf('Calculating splining matrix...\n')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Open mesh file - contains POS and index1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Project head vertices onto unit sphere
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
spherePOS = newPOS-ones(size(newPOS,1),1)*HeadCenter; % recenter
nPOSnorm = sqrt(sum(spherePOS.^2,2));
spherePOS(:,1) = spherePOS(:,1)./nPOSnorm;
spherePOS(:,2) = spherePOS(:,2)./nPOSnorm;
spherePOS(:,3) = spherePOS(:,3)./nPOSnorm;
x = spherePOS(:,1);
y = spherePOS(:,2);
z = spherePOS(:,3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate new electrode positions on head
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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); % closest 3 points
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 = newcoords;
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot mesh and electrodes only
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~strcmpi(g.plotmeshonly, 'off')
if strcmpi(g.plotmeshonly, 'sphere')
newElect(:,1) = Xe;
newElect(:,2) = Ye;
newElect(:,3) = Ze;
POS(index1,:) = spherePOS; HeadCenter = [ 0 0 0 ];
end
plotmesh(TRI1, POS, NORM);
plotelecopt.labelflag = 0;
plotelec(newElect, ElectrodeNames, HeadCenter, plotelecopt);
rotate3d;
return;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate g(x) for sphere mesh vertices
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Computing %d vertices. Should take a while (see wait bar)\n',...
length(x))
fprintf(' but does not have to be done again for this montage...\n');
icadefs;
gx = fastcalcgx(x,y,z,Xe,Ye,Ze);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Save spline file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
comment = g.comment;
headplot_version = 2;
transform = g.transform;
try, save(spline_file, '-V6', '-mat', 'Xe', 'Ye', 'Ze', 'G', 'gx', 'newElect', ...
'ElectrodeNames', 'indices', 'comment', 'headplot_version', 'transform');
catch,
try, save(spline_file, '-mat', 'Xe', 'Ye', 'Ze', 'G', 'gx', 'newElect', ...
'ElectrodeNames', 'indices', 'comment', 'headplot_version', 'transform');
catch, error('headplot: save spline file error, out of space or file permission problem');
end
end
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',...
'See also >> headplot cartesian\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',...
'See also >> headplot example\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] , [];
'lighting' , 'string', { 'on', 'off' } , 'on';
'verbose' 'string' { 'on','off' } 'on';
'maplimits' { 'string','real' } [] 'absmax';
'title' 'string' [] '';
'lights' 'real' [] DEFAULT_LIGHTS;
'view' { 'string','real' } [] [143 18];
'colormap' 'real' [] feval(DEFAULT_COLORMAP, 64);
'transform' 'real' [] [];
'meshfile' {'string','struct' } [] DEFAULT_MESH;
'electrodes' 'string' { 'on','off' } 'on';
'electrodes3d' 'string' { 'on','off' } 'off';
'eleccolor' 'cell' { } {};
'material' 'string' [] 'dull';
'orilocs' { 'string','struct' } [] '';
'labels' 'integer' [0 1 2] 0 }, 'headplot');
if ischar(g)
error(g);
end
plotelecopt.electrodes3d = g.electrodes3d;
if ~isempty(g.eleccolor) && length(g.eleccolor) ~= length(values)
error('The number of color must be the same as the number of channels to plot');
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');
if exist('indices')
try,
values = values(indices);
catch, error('problem of index or electrode number with splinefile'); end
end
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) = [ eloc_file(indices).X ]'; % attention inversion before
newElect(:,2) = [ eloc_file(indices).Y ]';
newElect(:,3) = [ eloc_file(indices).Z ]';
% optional transformation
% -----------------------
if ~isempty(g.transform)
transmat = traditionaldipfit( g.transform ); % arno
newElect = transmat*[ newElect ones(size(newElect,1),1)]';
newElect = newElect(1:3,:)';
end
end
% --------------
% load mesh file
% --------------
[newPOS, POS, TRI1, TRI2, NORM, index1, center] = getMeshData(g.meshfile);
%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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]);
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]);
else
ColorbarHandle = cbar(g.cbar,3,[amin amax]);
end
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)
try
lighting phong % all this gives a matte reflectance
catch
lighting gouraud % Octave
end
material(g.material);
end
%%%%%%%%%%%%%%%%%%%%%%%%%
% Set viewpoint
%%%%%%%%%%%%%%%%%%%%%%%%%
if ischar(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 ~ischar(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')
plotelecopt.labelflag = g.labels;
plotelecopt.eleccolor = g.eleccolor;
plotelec(newElect, ElectrodeNames, HeadCenter, plotelecopt);
else
fprintf('Variable newElect not read from spline file.\n');
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Turn on rotate3d, allowing rotation of the plot using the mouse
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ismatlab
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
else
rotate3d
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);
%%%%%%%%%%%%%%%%%%%
function gx = fastcalcgx(x,y,z,Xe,Ye,Ze)
onemat = ones(length(x),length(Xe));
EI = onemat - sqrt((repmat(x,1,length(Xe)) - repmat(Xe',length(x),1)).^2 +...
(repmat(y,1,length(Xe)) - repmat(Ye',length(x),1)).^2 +...
(repmat(z,1,length(Xe)) - repmat(Ze',length(x),1)).^2);
%
gx = zeros(length(x),length(Xe));
m = 4;
icadefs;
if ismatlab
hwb = waitbar(0,'Computing spline file (only done once)...', 'color', BACKEEGLABCOLOR);
hwbend = 7;
end
for n = 1:7
L = legendre(n,EI);
gx = gx + ((2*n+1)/(n^m*(n+1)^m))*squeeze(L(1,:,:));
if ismatlab
waitbar(n/hwbend,hwb);
end
end
gx = gx/(4*pi);
if ismatlab
close(hwb);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
% %%%%%%%%%%%%%%%
% plot electrodes
% %%%%%%%%%%%%%%%
function plotelec(newElect, ElectrodeNames, HeadCenter, opt);
newNames = newElect*opt.NamesDFac; % Calculate electrode label positions
for i = 1:size(newElect,1)
if newElect(i,:) ~= [0 0 0] % plot radial lines to electrode sites
if strcmpi(opt.electrodes3d, 'off')
line([newElect(i,1) HeadCenter(1)],[newElect(i,2) HeadCenter(2)],...
[newElect(i,3) HeadCenter(3)],'color',opt.MarkerColor,'linewidth',1);
end
if opt.labelflag == 1 % plot electrode numbers
t=text(newNames(i,1),newNames(i,2),newNames(i,3),int2str(i));
set(t,'Color',opt.NamesColor,'FontSize',opt.NamesSize,'FontWeight','bold',...
'HorizontalAlignment','center');
elseif opt.labelflag == 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',opt.NamesColor,'FontSize',opt.NamesSize,'FontWeight','bold',...
'HorizontalAlignment','center');
else
fprintf('Variable ElectrodeNames not read from spline file.\n');
end
else % plot electrode markers
if strcmpi(opt.electrodes3d, 'off')
if isempty(opt.eleccolor)
line(newElect(i,1),newElect(i,2),newElect(i,3),'marker',...
'.','markersize',20,'color',opt.MarkerColor,'linestyle','none');
else
if isempty(opt.eleccolor{i})
line(newElect(i,1),newElect(i,2),newElect(i,3),'marker',...
'.','markersize',12,'color',opt.MarkerColor,'linestyle','none');
else
line(newElect(i,1),newElect(i,2),newElect(i,3),'marker',...
'.','markersize',35,'color',opt.MarkerColor,'linestyle','none');
line(newElect(i,1),newElect(i,2),newElect(i,3),'marker',...
'.','markersize',28,'color',[1 1 1],'linestyle','none');
line(newElect(i,1),newElect(i,2),newElect(i,3),'marker',...
'.','markersize',15,'color',opt.eleccolor{i},'linestyle','none');
end
end
else
[xc yc zc] = cylinder( 2, 10);
[xs ys zs] = sphere(10);
xc = [ xc; -xs(7:11,:)*2 ];
yc = [ yc; -ys(7:11,:)*2 ];
zc = [ zc; zs(7:11,:)*0.2+1 ];
hold on;
cylinderSize = 3;
colorarray = repmat(reshape(opt.MarkerColor, 1,1,3), [size(zc,1) size(zc,2) 1]);
handles = surf(xc*cylinderSize, yc*cylinderSize, zc*cylinderSize, colorarray, 'edgecolor', 'none', ...
'backfacelighting', 'lit', 'facecolor', 'interp', 'facelighting', ...
'phong', 'ambientstrength', 0.3);
cylnderHeight = 1;
if newElect(i,3) < 10, addZ = -30; else addZ = 0; end
if newElect(i,3) < -20, addZ = -60; else addZ = 0; end
xx = newElect(i,1) - ( newElect(i,1)-HeadCenter(1) ) * 0.01 * cylnderHeight;
xxo1 = newElect(i,1) + ( newElect(i,1)-HeadCenter(1) ) * 0.01 * cylnderHeight;
yy = newElect(i,2) - ( newElect(i,2)-HeadCenter(2) ) * 0.01 * cylnderHeight;
yyo1 = newElect(i,2) + ( newElect(i,2)-HeadCenter(2) ) * 0.01 * cylnderHeight;
zz = newElect(i,3) - ( newElect(i,3)-HeadCenter(3)-addZ ) * 0.01 * cylnderHeight;
zzo1 = newElect(i,3) + ( newElect(i,3)-HeadCenter(3)-addZ ) * 0.01 * cylnderHeight;
[xc yc zc] = adjustcylinder2( handles, [xx yy zz], [xxo1 yyo1 zzo1] );
end
end
end
end
% get mesh information
% --------------------
function [newPOS POS TRI1 TRI2 NORM index1 center] = getMeshData(meshfile)
%#function mheadnew.mat
meshpath = fileparts(which('mheadnew.mat'));
if isdeployed
addpath( meshpath );
end
if ~isstruct(meshfile)
if ~exist(meshfile, 'file')
if isdeployed
meshfile = fullfile( meshpath, meshfile);
if ~exist(meshfile, 'file')
error('headplot(): deployed mesh file "%s" not found\n', meshfile);
end
else
error('headplot(): mesh file "%s" not found\n', meshfile);
end
end
fprintf('Loaded mesh file %s\n',meshfile);
try
meshfile = load(meshfile,'-mat');
catch
meshfile = [];
meshfile.POS = load('mheadnewpos.txt', '-ascii');
meshfile.TRI1 = load('mheadnewtri1.txt', '-ascii'); % upper head
%try, TRI2 = load('mheadnewtri2.txt', '-ascii'); catch, end; % lower head
%index1 = load('mheadnewindex1.txt', '-ascii');
meshfile.center = load('mheadnewcenter.txt', '-ascii');
end
end
if isfield(meshfile, 'vol')
if isfield(meshfile.vol, 'r')
[X Y Z] = sphere(50);
POS = { X*max(meshfile.vol.r) Y*max(meshfile.vol.r) Z*max(meshfile.vol.r) };
TRI1 = [];
else
POS = meshfile.vol.bnd(1).pnt;
TRI1 = meshfile.vol.bnd(1).tri;
end
elseif isfield(meshfile, 'bnd')
POS = meshfile.bnd(1).pnt;
TRI1 = meshfile.bnd(1).tri;
elseif isfield(meshfile, 'TRI1')
POS = meshfile.POS;
TRI1 = meshfile.TRI1;
try TRI2 = meshfile.TRI2; catch, end % NEW
try center = meshfile.center; catch, end % NEW
elseif isfield(meshfile, 'vertices')
POS = meshfile.vertices;
TRI1 = meshfile.faces;
else
error('Unknown Matlab mesh file');
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
newPOS = POS(index1,:);
Datasets
Models
