eeg_plot_surf_contours

 eeg_plot_surf_contours - plot 3D contours on triangulated surface
 
 Usage: [p] = eeg_contours_3d(p,[mode])
 
 p is the eeg_toolbox struct (see eeg_toolbox_defaults)

 mode is:

 'rb' for +ve red, -ve blue
 'bw' for +ve solid black, -ve dashed black (default)
 
 In development!

0001 function [ p ] = eeg_plot_surf_contours(p,mode)
0002 
0003 % eeg_plot_surf_contours - plot 3D contours on triangulated surface
0004 %
0005 % Usage: [p] = eeg_contours_3d(p,[mode])
0006 %
0007 % p is the eeg_toolbox struct (see eeg_toolbox_defaults)
0008 %
0009 % mode is:
0010 %
0011 % 'rb' for +ve red, -ve blue
0012 % 'bw' for +ve solid black, -ve dashed black (default)
0013 %
0014 % In development!
0015 %
0016 
0017 % $Revision: 1.1 $ $Date: 2004/11/12 01:32:33 $
0018 
0019 % Licence:  Gnu GPL, no express or implied warranties
0020 % History:  04/03 Darren.Weber_at_radiology.ucsf.edu
0021 %           obtained permission from Robert Oostenveld to
0022 %           distribute his code below under the GPL
0023 %
0024 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0025 
0026 
0027 fprintf('\nEEG_PLOT_SURF_CONTOURS (%s)...\n',['$Revision: 1.1 $']); tic;
0028 
0029 % confirm that a patch surface is available
0030 if isfield(p,'mesh'),
0031     if isfield(p.mesh,'data'),
0032         if isfield(p.mesh.data,'timeseries'),
0033             if isempty(p.mesh.data.Cdata{p.mesh.current}),
0034                 msg = sprintf('...p.mesh.data.Cdata{%d} is empty\n',p.mesh.current);
0035                 error(msg);
0036             end
0037         end
0038     end
0039 else
0040     error('...p.mesh.data is empty - load mesh first\n');
0041 end
0042 
0043 if ~exist('mode','var'), mode = 'bw'; end
0044 if isempty(mode),        mode = 'bw'; end
0045 
0046 
0047 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0048 % Determine/validate the min, max surface color data range
0049 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0050 
0051 % Select mesh timeseries values at samplePoint
0052 switch p.mesh.data.meshtype{p.mesh.current},
0053 case {'scalp','elec'},
0054     samplePoint = p.volt.samplePoint;
0055 otherwise
0056     samplePoint = p.mesh.samplePoint;
0057 end
0058 
0059 % get vertices and faces of mesh
0060 vertices = p.mesh.data.vertices{p.mesh.current};
0061 faces = p.mesh.data.faces{p.mesh.current};
0062 
0063 % get number of vertices
0064 nvert = size(p.mesh.data.vertices{p.mesh.current},1);
0065 % Assume more vertices than time points
0066 [s1,s2] = size(p.mesh.data.Cdata{p.mesh.current});
0067 if isequal(nvert,s1), % vertices in rows, timepoints in columns
0068     meshCdata = p.mesh.data.Cdata{p.mesh.current}(:,samplePoint);
0069 else,                 % vertices in columns, timeseries in rows
0070     meshCdata = p.mesh.data.Cdata{p.mesh.current}(samplePoint,:)';
0071 end
0072 
0073 switch p.rangeMethod
0074 case 'minmaxall', % Min/Max,all points
0075     fprintf('...estimating color data range, min/max all time points.\n');
0076     p.maximumIntensity = max(max(p.mesh.data.Cdata{p.mesh.current}));
0077     p.minimumIntensity = min(min(p.mesh.data.Cdata{p.mesh.current}));
0078 case 'minmaxone', % Min/Max, single point
0079     fprintf('...estimating color data range, min/max single time point.\n');
0080     % get number of vertices
0081     p.maximumIntensity = max(max(meshCdata));
0082     p.minimumIntensity = min(min(meshCdata));
0083 case 'minmaxabs', % Min/Max, Absolute
0084     fprintf('...estimating color data range, abs min/max single time point.\n');
0085     absmax = max(max(abs(meshCdata)));
0086     p.maximumIntensity =  absmax;
0087     p.minimumIntensity = -absmax;
0088 otherwise
0089     % check that specified intensity range is defined
0090     fprintf('...checking predefined color data range.\n');
0091     if isempty(p.maximumIntensity),
0092         fprintf('...estimating color data range, min/max single time point.\n');
0093         p.maximumIntensity = max(max(meshCdata)); end
0094     if isempty(p.minimumIntensity),
0095         fprintf('...estimating color data range, min/max single time point.\n');
0096         p.minimumIntensity = min(min(meshCdata)); end
0097 end
0098 
0099 
0100 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0101 % Determine the contour values
0102 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0103 
0104 switch p.contour.stepMethod
0105 case 0  % StepSize method
0106     fprintf('...using contour step size: %6.2f\n',p.contour.stepSize);
0107     p.contour.levels = eeg_contour_levels( p.contour.stepSize, meshCdata );
0108     p.contour.Nsteps = length(p.contour.levels);
0109 otherwise % Number of steps method
0110     fprintf('...using number of contours: %d\n',p.contour.Nsteps);
0111     p.contour.stepSize = abs(p.maximumIntensity - p.minimumIntensity) / p.contour.Nsteps;
0112     p.contour.levels = eeg_contour_levels( p.contour.stepSize, meshCdata );
0113 end
0114 
0115 
0116 % Exit if there are no contour levels
0117 if isempty(p.contour.levels),
0118     fprintf('...no contours in this data range.\n');
0119     % remove any current contours
0120     if isfield(p.contour,'patches'),
0121         if ~isempty(p.contour.patches),
0122             handleIndex = find(ishandle(p.contour.patches));
0123             delete(p.contour.patches(handleIndex));
0124         end
0125     end
0126     p.contour.patches = [];
0127     return
0128 end
0129 
0130 
0131 
0132 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0133 % This code below adapted from triplot.m (c) 2001 Robert Oostenveld ...
0134 % This code here is licenced under the GPL, as of April 8th 2003.
0135 
0136 fprintf('...calculating surface contours.\n');
0137 
0138 
0139 % create an Nx3 matrix of color data values at each face vertex
0140 FaceCdata = meshCdata(faces);
0141 % find the min/max vertex values of each face
0142 FaceCdata_min = min(FaceCdata, [], 2);
0143 FaceCdata_max = max(FaceCdata, [], 2);
0144 
0145 for contourIndex=1:length(p.contour.levels),
0146     
0147     % Find all the faces containing any vertices within the contour level
0148     contourLevel = p.contour.levels(contourIndex);
0149     use = contourLevel>=FaceCdata_min & contourLevel<=FaceCdata_max;
0150     faceIndices = find(use)';
0151     
0152     intersect1 = [];
0153     intersect2 = [];
0154     
0155     counter = 0;
0156     for faceIndex=faceIndices,
0157         
0158         xyz  = vertices(faces(faceIndex,:), :); % 3x3, rows = vertexIJK, columns = XYZ
0159         v = FaceCdata(faceIndex,:);             % 1x3, Cdata at rows of xyz
0160         
0161         % find line through face that passes through contourLevel ?
0162         
0163         la(1) = (contourLevel-v(1)) / (v(2)-v(1));    % abscissa between vertex 1 and 2
0164         la(2) = (contourLevel-v(2)) / (v(3)-v(2));    % abscissa between vertex 2 and 3
0165         la(3) = (contourLevel-v(3)) / (v(1)-v(3));    % abscissa between vertex 1 and 2
0166         
0167         abc(1,:) = xyz(1,:) + la(1) * (xyz(2,:) - xyz(1,:));
0168         abc(2,:) = xyz(2,:) + la(2) * (xyz(3,:) - xyz(2,:));
0169         abc(3,:) = xyz(3,:) + la(3) * (xyz(1,:) - xyz(3,:));
0170         
0171         counter = counter + 1;
0172         select  = find(la>=0 & la<=1);
0173         
0174         intersect1(counter, :) = abc(select(1),:);
0175         intersect2(counter, :) = abc(select(2),:);
0176     end
0177     
0178     % remember the details for external reference
0179     contour(contourIndex).level = contourLevel;
0180     contour(contourIndex).n     = counter;
0181     contour(contourIndex).intersect1 = intersect1;
0182     contour(contourIndex).intersect2 = intersect2;
0183 end
0184 
0185 
0186 % collect all different contourlevels
0187 intersect1 = [];
0188 intersect2 = [];
0189 contourLevel = [];
0190 for contourIndex=1:length(p.contour.levels)
0191     intersect1 = [intersect1; contour(contourIndex).intersect1];
0192     intersect2 = [intersect2; contour(contourIndex).intersect2];
0193     contourLevel = [contourLevel; ones(contour(contourIndex).n,1) * p.contour.levels(contourIndex)];
0194 end
0195 
0196 X = [intersect1(:,1) intersect2(:,1)]';
0197 Y = [intersect1(:,2) intersect2(:,2)]';
0198 C = [contourLevel(:) contourLevel(:)]';
0199 
0200 if size(vertices,2)>2
0201     % 3D contours
0202     Z = [intersect1(:,3) intersect2(:,3)]';
0203 else
0204     % 2D contours
0205     Z = zeros(2, length(contourLevel));
0206 end
0207 
0208 
0209 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0210 % Plot the surface contours
0211 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0212 
0213 % remove any current contours
0214 if isfield(p.contour,'patches'),
0215     if ~isempty(p.contour.patches),
0216         handleIndex = find(ishandle(p.contour.patches));
0217         delete(p.contour.patches(handleIndex));
0218     end
0219 end
0220 
0221 p.contour.patches = [];
0222 
0223 fprintf('...plotting surface contours.\n');
0224 
0225 for i=1:length(contourLevel)
0226     
0227     % make colour contours (default)
0228     h = patch('Tag','CONTOUR','XData', X(:,i), 'Ydata', Y(:,i), ...
0229         'ZData', Z(:,i), 'CData', C(:,i), ...
0230         'facecolor','none','edgecolor','flat', ...
0231         'linestyle', '-', 'linewidth', 1, ...
0232         'userdata',contourLevel(i));
0233     
0234     switch mode,
0235     case 'rb',
0236         % make red-blue contours
0237         if     contourLevel(i)>0,   edgecolor = 'red';
0238         elseif contourLevel(i)<0,   edgecolor = 'blue';
0239         else                        edgecolor = 'black';
0240         end
0241     otherwise
0242         % make black-white contours
0243         if     contourLevel(i)>0,   edgecolor = [.4 .4 .4];
0244         elseif contourLevel(i)<0,   edgecolor = [.2 .2 .2];
0245         else                        edgecolor = [.5 .5 .5];
0246         end
0247     end
0248     % apply contour styles
0249     if     contourLevel(i)>0,   set(h,'edgecolor',edgecolor,'linestyle', '-');
0250     elseif contourLevel(i)<0,   set(h,'edgecolor',edgecolor,'linestyle', ':','linewidth', 1);
0251     else                        set(h,'edgecolor',edgecolor,'linewidth', 1.5);
0252     end
0253     
0254     p.contour.patches = [p.contour.patches; h];
0255 end
0256 
0257 t = toc;
0258 fprintf('...done (%6.2f sec).\n\n',t);
0259 
0260 return