eeg_contours_engine

 eeg_contours_engine - Generates contour/topographic maps
 
 Usage: [p] = eeg_contours_engine(p)

 p is the eeg_toolbox struct (see eeg_toolbox_defaults)

 Given:
   data file of electrode positions (see elec_load for help)
   data file of electrode voltage values (see eeg_ascii_load for help)
   other parameters for selection and plot control

   Parameters can be passed as a structure (p) with the function call or
   passed in the file 'eeg_contours_default.txt'.  A standard set of
   parameters can be created using 'eeg_contours_create_defaults'. See
   that function for more information about the parameters created.

 Output:
    various topographic maps and 2D or 3D contour plots

0001 function [p] = eeg_contours_engine(p)
0002 
0003 % eeg_contours_engine - Generates contour/topographic maps
0004 %
0005 % Usage: [p] = eeg_contours_engine(p)
0006 %
0007 % p is the eeg_toolbox struct (see eeg_toolbox_defaults)
0008 %
0009 % Given:
0010 %   data file of electrode positions (see elec_load for help)
0011 %   data file of electrode voltage values (see eeg_ascii_load for help)
0012 %   other parameters for selection and plot control
0013 %
0014 %   Parameters can be passed as a structure (p) with the function call or
0015 %   passed in the file 'eeg_contours_default.txt'.  A standard set of
0016 %   parameters can be created using 'eeg_contours_create_defaults'. See
0017 %   that function for more information about the parameters created.
0018 %
0019 % Output:
0020 %    various topographic maps and 2D or 3D contour plots
0021 %
0022 
0023 % $Revision: 1.2 $ $Date: 2005/07/13 06:02:11 $
0024 
0025 % Licence:  Gnu GPL, no express or implied warranties
0026 % History:  08/99 Chris Harvey,    Created, combining work by Darren Weber, Murk Bottomer
0027 %           10/99 CRH    Accepts input parameters, provides some validation
0028 %           07/01 Darren.Weber_at_radiology.ucsf.edu
0029 %               - modified parameter handling, now using 'p' structure
0030 %               - now calls separate functions to create/read/write defaults
0031 %               - modified most aspects of code to adapt to development of other
0032 %                 functions called from this function
0033 % Develop:  08/01 Darren.Weber_at_radiology.ucsf.edu
0034 %               - Needs scripting facility to wind through time points in
0035 %                 a range of timepoints, outputing contour plot(s) for each in
0036 %                 a given image format (EPS, GIF, TIF, etc.) or a movie of
0037 %                 the timeseries.
0038 %           05/02 Darren.Weber_at_radiology.ucsf.edu
0039 %               - now integrated with mesh interpolation and calls
0040 %                 animation interface, with image saving options
0041 %
0042 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0043 
0044 if ~exist('p','var'),
0045   fprintf('Setting default parameters.\n');
0046  [p] = eeg_toolbox_defaults('create');
0047 end
0048 
0049 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0050 % load electrode co-ordinates
0051 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0052 
0053 if isempty(p.elec.data),[p] = elec_open(p); end
0054 
0055 %elecLabels = p.elec.data.label;
0056 X = p.elec.data.x;
0057 Y = p.elec.data.y;
0058 Z = p.elec.data.z;
0059 
0060 Xrad = p.elec.data.R(1);
0061 Yrad = p.elec.data.R(2);
0062 Zrad = p.elec.data.R(3);
0063 
0064 elecNumElectrodes = length(p.elec.data.x);
0065 
0066 % Estimate sphere or ellipse that best fits electrode co-ordinates
0067 switch p.elec.shape
0068   case 'ellipse'
0069     Xp = p.elec.data.Xel;
0070     Yp = p.elec.data.Yel;
0071     Zp = p.elec.data.Zel;
0072   case 'sphere'
0073     Xp = p.elec.data.Xsp;
0074     Yp = p.elec.data.Ysp;
0075     Zp = p.elec.data.Zsp;
0076   otherwise
0077     Xp = p.elec.data.x;
0078     Yp = p.elec.data.y;
0079     Zp = p.elec.data.z;           
0080 end
0081 
0082 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0083 % Read voltage data
0084 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0085 
0086 if isempty(p.volt.data),[p] = eeg_open(p); end
0087 
0088 [voltNumTimePoints, voltNumElectrodes] = size(p.volt.data);
0089 
0090 % Validate voltage matrix orientation against electrode file
0091 
0092 if (voltNumElectrodes ~= elecNumElectrodes)
0093   fprintf('\nWarning: # electrodes, data file %d, electrode file %d\n', voltNumElectrodes, elecNumElectrodes);
0094   if (voltNumTimePoints == elecNumElectrodes)
0095     % Assume data file needs rotation
0096     fprintf('Continuing with data file rotated\n');
0097     p.volt.data = p.volt.data';   [voltNumTimePoints, voltNumElectrodes] = size(p.volt.data);
0098   else
0099     fprintf('\nError: Cannot reconcile electrodes with voltage data.\n');
0100     return
0101   end
0102 end
0103 
0104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0105 % Select a timePoint
0106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0107 if or(isequal(p.clickTimePoint, 1),isempty(p.volt.samplePoint))
0108   figure('NumberTitle','off','Name','MouseClick the TimePoint')
0109   plot(p.volt.data); 
0110   [Xt,Yt] = ginput(1);
0111   p.volt.samplePoint = round(Xt);
0112   close
0113 end
0114 if (p.volt.samplePoint > voltNumTimePoints)
0115   msg = sprintf(' Selected timepoint exceeds data points of %d\n\n', voltNumTimePoints);
0116   warning(msg);
0117   p.volt.samplePoint = voltNumTimePoints;
0118 end
0119 if ~isempty(p.volt.timeArray),
0120   p.volt.sampleTime = p.volt.timeArray(p.volt.samplePoint);
0121 end
0122 
0123 V = p.volt.data(p.volt.samplePoint,:)';
0124 
0125 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0126 % Determine/validate the min, max range and the step size
0127 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0128 
0129 switch p.rangeMethod
0130   case 'minmaxall', % Min/Max,all points
0131     p.maximumIntensity = max(max(p.volt.data));
0132     p.minimumIntensity = min(min(p.volt.data));
0133   case 'minmaxone', % Min/Max, single point
0134     p.maximumIntensity = max(max(V));
0135     p.minimumIntensity = min(min(V));
0136   case 'minmaxabs', % Min/Max, Absolute
0137     absmax = max(max(abs(V)));
0138     p.maximumIntensity =  absmax;
0139     p.minimumIntensity = -absmax;
0140   otherwise
0141     % check that specified intensity range is defined
0142     if isempty(p.maximumIntensity),
0143       p.maximumIntensity = max(max(V)); end
0144     if isempty(p.minimumIntensity),
0145       p.minimumIntensity = min(min(V)); end
0146 end
0147 
0148 switch p.contour.stepMethod
0149   case 0  % StepSize method
0150     p.contour.levels = eeg_contour_levels( p.contour.stepSize, V );
0151     p.contour.Nsteps = length(p.contour.levels);
0152   otherwise % Number of steps method
0153     p.contour.stepSize = abs(p.maximumIntensity - p.minimumIntensity) / p.contour.Nsteps;
0154     p.contour.levels = eeg_contour_levels( p.contour.stepSize, V );
0155 end
0156 
0157 
0158 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0159 % Draw the plots
0160 %    X,Y,Z            = actual electrode placements
0161 %     Xel,Yel,Zel     = electrode position on the ellipsoid surface
0162 %    gridSize        = number of grid lines to use on each plane
0163 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0164 
0165 
0166 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0167 % Interpolate surface
0168 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0169 
0170 % Determine interpolation gridding
0171 m = max(abs([Xp;Yp]));
0172 m = ceil(m);
0173 if (p.grid.method == 2), xg = -m:p.grid.res:m;             % grid resolution method
0174 else,                   xg = linspace(-m,m,p.grid.size);  % grid size method
0175 end
0176 yg = xg';
0177 [Xi,Yi,Zi] = griddata(Xp,Yp,Zp,xg,yg,p.interpMethod);    % Surface interpolation
0178 Vi         = griddata(Xp,Yp,V ,xg,yg,p.interpMethod);    % Voltage interpolation
0179 
0180 % Calculate voltage at refined electrode vertices
0181 %if ~isempty(p.elec.data.Vsp),
0182 %    p.elec.data.Psp = p.elec.data.Isp * V;
0183 %end
0184 
0185 
0186 % Calculate voltage at scalp mesh vertices
0187 if p.mesh.plotSurf,
0188   
0189   % get the scalp mesh
0190   [p.mesh.current,meshExists] = mesh_check(p,'scalp');
0191   
0192   if isempty(p.mesh.current),
0193     fprintf('...no scalp mesh, using electrode surface.\n');
0194     p.mesh.plotSurf = 0;
0195     p.elec.plotSurf = 1;
0196   else
0197     % Assume that p.mesh.data has this field, as it should
0198     % be initialised by mesh_open, but check anyway
0199     if ~isfield(p.mesh.data,'Cdata'),
0200      [p] = mesh_scalp_interp(p);
0201     end
0202     % Confirm that p.mesh.data.timeseries{p.mesh.current} has data
0203     if length(p.mesh.data.Cdata) < p.mesh.current,
0204      [p] = mesh_scalp_interp(p);
0205     end
0206     if isempty(p.mesh.data.Cdata{p.mesh.current}),
0207      [p] = mesh_scalp_interp(p);
0208     end
0209     
0210     % Now we should have a mesh timeseries
0211     
0212     % Verify the timeseries data against the elec timeseries;
0213     % if generated with mesh_scalp_interp, they should be equal
0214     TMP = p.mesh.data.Cdata{p.mesh.current};
0215     Nelec = size(p.volt.data,2);
0216     if ~isequal(p.volt.data',TMP(1:Nelec,:)),
0217       % Something has changed, either a new electrode
0218       % set or scalp mesh, so run the interpolation
0219       msg = sprintf('scalp time series not consistent with electrode time series - recalculating.\n');
0220       warning(msg);
0221       p.mesh.data.Cdata{p.mesh.current} = [];
0222      [p] = mesh_scalp_interp(p);
0223       TMP = p.mesh.data.Cdata{p.mesh.current};
0224       if ~isequal(p.volt.data',TMP(1:Nelec,:)),
0225         msg = sprintf('scalp timeseries not consistent with electrode timeseries.\n');
0226         error(msg);
0227       end
0228     end
0229     clear TMP;
0230   end
0231 end
0232 
0233 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0234 % surface plot
0235 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0236 
0237 if or(p.elec.plotSurf,p.mesh.plotSurf),
0238   
0239   if ~isequal(p.volt.sampleTime,0),
0240     name = sprintf('Surface: %s @ %8.2f msec',p.volt.file, p.volt.sampleTime);
0241   else
0242     name = sprintf('Surface: %s',p.volt.file);
0243   end
0244   
0245   if p.elec.plotSurf,
0246     
0247     [p.mesh.current,meshExists] = mesh_check(p,'elec');
0248     
0249     if ~meshExists,
0250       p.mesh.data.meshtype{p.mesh.current}    = 'elec';
0251       p.mesh.data.vertices{p.mesh.current}    = [Xp,Yp,Zp];
0252       p.mesh.data.faces{p.mesh.current}       = convhulln([Xp,Yp,Zp]);
0253       p.mesh.data.lapmat{p.mesh.current}      = [];
0254       p.mesh.data.lapint{p.mesh.current}      = [];
0255       p.mesh.data.timeseries{p.mesh.current}  = p.volt.timeArray(:,1);
0256       p.mesh.data.Cdata{p.mesh.current}       = p.volt.data';
0257     end
0258     
0259     %p = eeg_plot_surf(p);
0260     %return
0261     
0262     H.gui = figure('NumberTitle','off','Name',name,...
0263       'PaperType','A4','PaperUnits','centimeters');
0264     
0265     set(gca,'Projection','perspective')
0266     set(gca,'DataAspectRatio',[1 1 1]);
0267     
0268     p.elec.patch = patch('vertices',[Xp,Yp,Zp],'faces',convhulln([Xp,Yp,Zp]),...
0269       'facevertexCdata',V,'facecolor','interp','EdgeColor','none',...
0270       'FaceLighting','phong');
0271     %patch('vertices',p.elec.data.Vsp,'faces',p.elec.data.Fsp,...
0272     %      'facevertexCdata',p.elec.data.Psp,'facecolor','interp','EdgeColor','none');
0273     
0274     set(gca,'Projection','perspective')
0275     set(gca,'DataAspectRatio',[1 1 1]);
0276     axis off tight vis3d
0277     
0278     lighting phong
0279     set(H.gui,'Renderer','zbuffer')
0280     
0281     % Create a light above the z axis
0282     H.light = light('style','infinite','position',[0 0 1000]);
0283     
0284     % For a near skin color
0285     %set(gca,'AmbientLightColor',[.9 .8 .7]);
0286     
0287     % MATERIAL([ka kd ks n sc]) sets the ambient/diffuse/specular strength,
0288     %    specular exponent and specular color reflectance of the objects.
0289     %p.reflect = material('dull');
0290     p.reflect{1} = 0.9;
0291     p.reflect{2} = 0.1;
0292     p.reflect{3} = 0.0;
0293     p.reflect{4} = 500;
0294     p.reflect{5} = 0;
0295     set(p.elec.patch,'AmbientStrength',p.reflect{1});
0296     set(p.elec.patch,'DiffuseStrength',p.reflect{2});
0297     set(p.elec.patch,'SpecularStrength',p.reflect{3});
0298     set(p.elec.patch,'SpecularExponent',p.reflect{4});
0299     set(p.elec.patch,'SpecularColorReflectance',p.reflect{5});
0300     
0301     if p.elec.plot,
0302       hold on, plot3(Xp,Yp,Zp,'k.');
0303     end
0304     
0305     set(gca,'Visible','off');
0306     colormap(p.colorMap.map);
0307     caxis([p.minimumIntensity p.maximumIntensity]);
0308     colorbar
0309     
0310     H.p = p;
0311     H.p.mesh.plotSurf = 0;
0312     set(H.gui,'userdata',H);
0313     
0314     SaveGraphics(H.gui,'_3Delec_',p);
0315     if isequal(exist('mouse_rotate'),2),
0316       mouse_rotate;
0317     else
0318       rotate3D;
0319     end
0320     if isequal(exist('gui_topo_animate'),2),
0321       gui_topo_animate('init',p);
0322     end
0323   end
0324   
0325   
0326   % Plot a mesh surface
0327   
0328   if p.mesh.plotSurf,
0329     
0330     [p.mesh.current,meshExists] = mesh_check(p,'scalp');
0331     
0332    [p] = eeg_plot_surf(p);
0333     
0334     return
0335   end
0336 end
0337 
0338 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0339 % Contour in 2-D
0340 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0341 
0342 %CONTOUR(Z,N) and CONTOUR(X,Y,Z,N) draw N contour lines,
0343 %overriding the automatic value.
0344 %CONTOUR(Z,V) and CONTOUR(X,Y,Z,V) draw LENGTH(V) contour lines
0345 %at the values specified in vector V.
0346 
0347 if isequal(p.contour.raw2D,1),
0348   
0349   if ~isequal(p.volt.sampleTime,0),
0350     name = sprintf('2-D Contour: %s @ %8.2f msec',p.volt.file, p.volt.sampleTime);
0351   else
0352     name = sprintf('2-D Contour: %s',p.volt.file);
0353   end
0354   fig = figure('NumberTitle','off','Name',name); colormap(p.colorMap.map);
0355   set(gca,'Projection','perspective')
0356   %set(gca,'Projection','orthographic')
0357   set(gca,'DataAspectRatio',[1 1 1]);
0358   
0359   %trisurf(delaunay(Xp,Yp),Xp,Yp,Zp,V,'EdgeColor','none','FaceColor','interp');
0360   %hold on, plot3(Xp,Yp,Zp,'.');
0361   %contour(Xi,Yi,Vi, p.contour.Nsteps) %, colorbar, hold on
0362   
0363   [c,ch,cf] = contourf(Xi,Yi,Vi,p.contour.levels);
0364   absmax = max(max(abs(Vi)));caxis([-absmax absmax]);
0365   %[C,CH,CF] = CONTOURF(...) also returns a column vector CH of handles
0366   %to PATCH objects and the contour matrix CF for the filled areas. The
0367   %UserData property of each object contains the height value for each
0368   %contour.
0369   
0370   % FaceColor is the contour patch color
0371   % EdgeColor is the contour line color
0372   for i=1:size(ch),
0373     height = get(ch(i),'UserData');
0374     
0375     if     (height > 0),
0376       set(ch(i),'LineStyle', '-'); % options are: -, --, :, -., none
0377       set(ch(i),'EdgeColor',[1 1 1]) % white contour lines, black = [0 0 0] (ie, [R G B])
0378     elseif (height < 0),
0379       set(ch(i),'LineStyle', ':'); % options are: -, --, :, -., none
0380     else
0381       set(ch(i),'LineStyle', '-'); % options are: -, --, :, -., none
0382       set(ch(i),'LineWidth',2);
0383     end
0384   end
0385   
0386   set(gca,'Visible','off'); %hold on; DrawHead(Xrad,'black'); hold off;
0387   colorbar
0388   
0389   SaveGraphics(fig,'_2D_',p);
0390   if isequal(exist('mouse_rotate'),2),
0391     mouse_rotate;
0392   else
0393     rotate3D;
0394   end
0395 end
0396 
0397 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0398 % Projected Contour in 2-D
0399 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0400 
0401 if isequal(p.contour.plot2D,1)
0402   
0403   % Project 3D electrode locations onto 2D plane.
0404   [X2p,Y2p] = elec_3d_2d(Xp,Yp,Zp,Zrad);
0405   m = max(abs([X2p;Y2p]));
0406   m = ceil(m);
0407   if (p.grid.method == 2)
0408     % Use grid resolution method
0409     X2g = -m:p.grid.res:m;
0410   else
0411     % Use grid size method
0412     X2g = linspace( -m, m, p.grid.size );
0413   end
0414   Y2g = X2g';
0415   [X2i Y2i V2i] = griddata(X2p,Y2p,V,X2g,Y2g,p.interpMethod);
0416   
0417   if ~isequal(p.volt.sampleTime,0),
0418     name = sprintf('Projected 2D Contours: %s @ %8.2f msec',p.volt.file, p.volt.sampleTime);
0419   else
0420     name = sprintf('Projected 2D Contours: %s',p.volt.file);
0421   end
0422   fig = figure('NumberTitle','off','Name',name);
0423   set(gca,'Projection','perspective')
0424   %set(gca,'Projection','orthographic')
0425   set(gca,'DataAspectRatio',[1 1 1]);
0426   
0427   colormap(p.colorMap.map);
0428   contourf(X2i,Y2i,V2i, p.contour.Nsteps);
0429   absmax = max(max(abs(V2i)));caxis([-absmax absmax]);colorbar
0430   xlabel('X'); ylabel('Y');
0431   
0432   SaveGraphics(fig,'_proj2D_',p);
0433 end
0434 
0435 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0436 % Projected Contour in 3-D
0437 % Draw the contour lines,
0438 %   Get the contours in 2-D.
0439 %   For each point in each contour line, know x,y
0440 %   Calculate z for the x,y from ellipse formula
0441 %   Plot3 point x,y,z
0442 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0443 
0444 if (p.contour.plot3D == 1)
0445   
0446   
0447   % Project 3-D electrode locations onto 2-D plane.
0448   [X2p,Y2p] = elec_3d_2d(Xp,Yp,Zp,Zrad);
0449   m = max(abs([X2p;Y2p]));
0450   m = ceil(m);
0451   if isequal(p.grid.method,2)
0452     % Use grid resolution method
0453     X2g = -m:p.grid.res:m;
0454   else
0455     % Use grid size method
0456     X2g = linspace( -m, m, p.grid.size );
0457   end
0458   Y2g = X2g';
0459   
0460   [X2i Y2i V3i] = griddata(Xp, Yp, V,X2g,Y2g,p.interpMethod);
0461   
0462   % Caculate the contour matrix
0463   C = contourc(X2g, Y2g, V3i, p.contour.Nsteps);
0464   
0465   if ~isequal(p.volt.sampleTime,0),
0466     name = sprintf('3D Contours: %s @ %8.2f msec',p.volt.file, p.volt.sampleTime);
0467   else
0468     name = sprintf('3D Contours: %s',p.volt.file);
0469   end
0470   
0471   %fig = figure('NumberTitle','off','Name',name);
0472   %set(gca,'Projection','perspective')
0473   %set(gca,'DataAspectRatio',[1 1 1]);
0474   %[C,Cpatch] = contourf(X2g, Y2g, V3i, p.contour.Nsteps);
0475   %get(Cpatch(1))
0476   %P.cdata    = get(Cpatch,'CData');
0477   %P.vertices = get(Cpatch,'Vertices');
0478   %P.faces    = get(Cpatch,'Faces');
0479   %close(fig);
0480   %fig = figure('NumberTitle','off','Name',name); hold on
0481   %i = 0;
0482   %while i < (length(Cpatch) - 1),
0483   %%vertices = P.vertices{i};
0484   %i = i + 1;
0485   %vertices = [P.vertices{i}; P.vertices{i+1}];
0486   %Xc = vertices(:,1);
0487   %Yc = vertices(:,2);
0488   %[X3,Y3,Z3] = elec_2d_3d(Xc,Yc,Xrad,Yrad,Zrad);
0489   %tri = delaunay(X3,Y3);
0490   %Htri = trisurf(tri,X3,Y3,Z3,P.cdata{i},'EdgeColor','none','FaceColor','flat');
0491   
0492   %%vertices = [X3 Y3 Z3];
0493   %%Hpatch = patch('Vertices',vertices,'Faces',P.faces{i},'CData',P.cdata{i},'FaceColor','flat');
0494   %contourValues(i) = P.cdata{i};
0495   %end
0496   
0497   %view(0,90), axis tight, rotate3d
0498   %set(gca,'Visible','off');
0499   %colormap(p.colorMap.map);
0500   %absmax = max(abs(contourValues));
0501   %caxis([-absmax absmax]);colorbar;
0502   
0503   
0504   
0505   
0506   fig = figure('NumberTitle','off','Name',name); hold on
0507   set(gca,'Projection','perspective')
0508   %set(gca,'Projection','orthographic')
0509   set(gca,'DataAspectRatio',[1 1 1]);
0510   
0511   % When we view the contour lines, those on the other side of
0512   % the head can be seen, which makes it confusing.
0513   
0514   % The array C is 2-row, x on 1st, y on 2nd.
0515   % Each contour set has a preceeding column, where
0516   % 1st is value of contour and 2nd is number of points.
0517   %
0518   % We need to get each contour separately, and calculate
0519   % Z for each X,Y.
0520   i = 1; j = 0;
0521   limit = size(C,2);
0522   while (i < limit),
0523     
0524     j = j + 1; contourValues(j) = C(1,i);
0525     
0526     numberPoints = C(2,i);
0527     
0528     endContourSet = i + numberPoints;
0529     
0530     Xc = [ C(1,i+1:endContourSet) ]';
0531     Yc = [ C(2,i+1:endContourSet) ]';
0532     
0533     %cdata = contourValues(j) + 0*Xc;  % Make cdata the same size as xdata
0534     
0535     % Given each X,Y we project back from 2D to 3D.
0536     [X3,Y3,Z3] = elec_2d_3d(Xc,Yc,Xrad,Yrad,Zrad);
0537     
0538     if (contourValues(j) < 0)
0539       % Values < 0 with dashed black line
0540       if isequal(p.colorMap.style,'Gray')
0541         line(X3,Y3,Z3,'LineStyle',':','LineWidth',0.5,'Color','black');
0542       else
0543         line(X3,Y3,Z3,'LineStyle',':','LineWidth',0.5,'Color','blue');
0544       end
0545     elseif (contourValues(j) > 0)
0546       % Values > 0 with solid black line
0547       if isequal(p.colorMap.style,'Gray')
0548         line(X3,Y3,Z3,'LineStyle','-','LineWidth',0.5,'Color','black');
0549       else
0550         line(X3,Y3,Z3,'LineStyle','-','LineWidth',0.5,'Color','red');
0551       end
0552     else
0553       % Values = 0 with thicker solid black line
0554       line(X3,Y3,Z3,'LineStyle','-','LineWidth',0.75,'Color','black');
0555     end
0556     i = endContourSet + 1;
0557   end
0558   
0559   view(0,90), axis tight, rotate3d
0560   set(gca,'Visible','off');
0561   %colormap(p.colorMap.map); H = colorbar;
0562   %absmax = max(abs(contourValues));
0563   
0564   SaveGraphics(fig,'_3D_',p);
0565 end
0566 
0567 return    
0568 
0569 
0570 
0571 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0572 function DrawHead(Hradius, Hcolor)
0573 
0574 switch Hcolor
0575   case 'black'
0576     HCOLOR = [0 0 0];
0577   case 'white'
0578     HCOLOR = [1 1 1];
0579   otherwise
0580     HCOLOR = [0 0 0];
0581 end
0582 HLINEWIDTH = 30;
0583 
0584 rmax = 0.95 * Hradius;
0585 
0586 % Plot Head
0587 
0588 theta1 = linspace(0,2*pi,50);       % 360 degree rotation
0589 thetad = (theta1(2) - theta1(1))/2; % provide staggered rotation
0590 theta = [theta1, theta1 + thetad];
0591 
0592 plot(cos(theta).*rmax,sin(theta).*rmax,'color',HCOLOR,'LineWidth',HLINEWIDTH);
0593 
0594 % Plot Nose
0595 width = rmax * 0.1;
0596 tip   = rmax * 1.075;
0597 base  = rmax * 1.005;
0598 
0599 plot([width;0;-width],[rmax;tip;rmax],'Color',HCOLOR,'LineWidth',HLINEWIDTH);
0600 return
0601 
0602 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0603 function SaveGraphics(F,type,p)
0604 if ~isequal(p.saveGraphics,'No Save Plots'),
0605   
0606   [path,file,ext] = fileparts(strcat(p.volt.path, filesep, p.volt.file));
0607   file = strcat(file, type, num2str(p.volt.samplePoint),'.',p.saveGraphics);
0608   file = fullfile(path,file);
0609   saveas(F,file);
0610 end
0611 return