eeg_contours_engine.m

Matlab下的EEG处理程序库

function [p] = eeg_contours_engine(p)
% EEG_CONTOURS_ENGINE - Generates contour/topographic maps
% 
% Useage:   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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% $Revision: 1.4 $ $Date: 2003/04/07 06:12:02 $

% Licence:  Gnu GPL, no express or implied warranties
% History:  08/99 Chris Harvey,	Created, combining work by Darren Weber, Murk Bottomer
%           10/99 CRH	Accepts input parameters, provides some validation
%           07/01 [e-mail]
%               - modified parameter handling, now using 'p' structure
%               - now calls separate functions to create/read/write defaults
%               - modified most aspects of code to adapt to development of other
%                 functions called from this function
% Develop:  08/01 [e-mail]
%               - Needs scripting facility to wind through time points in
%                 a range of timepoints, outputing contour plot(s) for each in
%                 a given image format (EPS, GIF, TIF, etc.) or a movie of
%                 the timeseries.
%           05/02 [e-mail]
%               - now integrated with mesh interpolation and calls
%                 animation interface, with image saving options
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('p','var'),
    fprintf('Setting default parameters.\n');
    p = eeg_toolbox_defaults('create');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load electrode co-ordinates 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if isempty(p.elec.data), p = elec_open(p); end
    
    elec = p.elec.data.label;
    X = p.elec.data.x;
    Y = p.elec.data.y;
    Z = p.elec.data.z;
    
    Xrad = p.elec.data.R(1);
    Yrad = p.elec.data.R(2);
    Zrad = p.elec.data.R(3);
    
    elecNumElectrodes = length(p.elec.data.x);
    
    % Estimate sphere or ellipse that best fits electrode co-ordinates
    switch p.elec.shape
       case 'ellipse'
           Xp = p.elec.data.Xel;
           Yp = p.elec.data.Yel;
           Zp = p.elec.data.Zel;
       case 'sphere'
           Xp = p.elec.data.Xsp;
           Yp = p.elec.data.Ysp;
           Zp = p.elec.data.Zsp;
       otherwise
           Xp = p.elec.data.x;
           Yp = p.elec.data.y;
           Zp = p.elec.data.z;           
    end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read voltage data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if isempty(p.volt.data), p = eeg_open(p); end
    
    [voltNumTimePoints, voltNumElectrodes] = size(p.volt.data);

    % Validate voltage matrix orientation against electrode file
    
    if (voltNumElectrodes ~= elecNumElectrodes)
    	fprintf('\nWarning: # electrodes, data file %d, electrode file %d\n', voltNumElectrodes, elecNumElectrodes);
    	if (voltNumTimePoints == elecNumElectrodes)
    		% Assume data file needs rotation
    		fprintf('Continuing with data file rotated\n');
    		p.volt.data = p.volt.data';   [voltNumTimePoints, voltNumElectrodes] = size(p.volt.data);
    	else
        	fprintf('\nError: Cannot reconcile electrodes with voltage data.\n');
    		return
    	end
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Select a timePoint
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if or(isequal(p.clickTimePoint, 1),isempty(p.volt.samplePoint))
    	figure('NumberTitle','off','Name','MouseClick the TimePoint')
    	plot(p.volt.data); 
    	[Xt,Yt] = ginput(1);
    	p.volt.samplePoint = round(Xt);
        close
    end
    if (p.volt.samplePoint > voltNumTimePoints)
    	msg = sprintf(' Selected timepoint exceeds data points of %d\n\n', voltNumTimePoints);
        warning(msg);
        p.volt.samplePoint = voltNumTimePoints;
    end
    if ~isempty(p.volt.timeArray),
        p.volt.sampleTime = p.volt.timeArray(p.volt.samplePoint);
    end
    
    V = p.volt.data(p.volt.samplePoint,:)';
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine/validate the min, max range and the step size
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    switch p.rangeMethod
    case 'minmaxall', % Min/Max,all points
    	p.maximumIntensity = max(max(p.volt.data));
        p.minimumIntensity = min(min(p.volt.data));
    case 'minmaxone', % Min/Max, single point
    	p.maximumIntensity = max(max(V));
        p.minimumIntensity = min(min(V));
    case 'minmaxabs', % Min/Max, Absolute
        absmax = max(max(abs(V)));
    	p.maximumIntensity =  absmax;
        p.minimumIntensity = -absmax;
    otherwise
        % check that specified intensity range is defined
        if isempty(p.maximumIntensity),
            p.maximumIntensity = max(max(V)); end
        if isempty(p.minimumIntensity),
            p.minimumIntensity = min(min(V)); end
    end
    
    switch p.contour.stepMethod
    case 0  % StepSize method
        p.contour.levels = eeg_contour_levels( p.contour.stepSize, V );
        p.contour.Nsteps = length(p.contour.levels);
    otherwise % Number of steps method
    	p.contour.stepSize = abs(p.maximumIntensity - p.minimumIntensity) / p.contour.Nsteps;
        p.contour.levels = eeg_contour_levels( p.contour.stepSize, V );
    end
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Draw the plots
%	X,Y,Z		    = actual electrode placements
% 	Xel,Yel,Zel 	= electrode position on the ellipsoid surface
%	gridSize	    = number of grid lines to use on each plane
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Interpolate surface
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Determine interpolation gridding
    m = max(abs([Xp;Yp]));
    m = ceil(m);
    if (p.grid.method == 2), xg = -m:p.grid.res:m;             % grid resolution method
    else,                   xg = linspace(-m,m,p.grid.size);  % grid size method
    end
    yg = xg';
    [Xi,Yi,Zi] = griddata(Xp,Yp,Zp,xg,yg,p.interpMethod);    % Surface interpolation
    Vi         = griddata(Xp,Yp,V ,xg,yg,p.interpMethod);    % Voltage interpolation
    
    % Calculate voltage at refined electrode vertices
    %if ~isempty(p.elec.data.Vsp),
    %    p.elec.data.Psp = p.elec.data.Isp * V;
    %end
    
    
    % Calculate voltage at scalp mesh vertices
    if p.mesh.plotSurf,
        
        % get the scalp mesh
        [p.mesh.current,meshExists] = mesh_check(p,'scalp');
        
        if isempty(p.mesh.current),
            fprintf('...no scalp mesh, using electrode surface.\n');
            p.mesh.plotSurf = 0;
            p.elec.plotSurf = 1;
        else
            % Assume that p.mesh.data has this field, as it should
            % be initialised by mesh_open, but check anyway
            if ~isfield(p.mesh.data,'Cdata'),
                p = mesh_scalp_interp(p);
            end
            % Confirm that p.mesh.data.timeseries{p.mesh.current} has data
            if length(p.mesh.data.Cdata) < p.mesh.current,
                p = mesh_scalp_interp(p);
            end
            if isempty(p.mesh.data.Cdata{p.mesh.current}),
                p = mesh_scalp_interp(p);
            end
            
            % Now we should have a mesh timeseries
            
            % Verify the timeseries data against the elec timeseries;
            % if generated with mesh_scalp_interp, they should be equal
            TMP = p.mesh.data.Cdata{p.mesh.current};
            Nelec = size(p.volt.data,2);
            if ~isequal(p.volt.data',TMP(1:Nelec,:)),
                % Something has changed, either a new electrode 
                % set or scalp mesh, so run the interpolation
                msg = sprintf('scalp time series not consistent with electrode time series - recalculating.\n');
                warning(msg);
                p.mesh.data.Cdata{p.mesh.current} = [];
                p = mesh_scalp_interp(p);
                TMP = p.mesh.data.Cdata{p.mesh.current};
                if ~isequal(p.volt.data',TMP(1:Nelec,:)),
                    msg = sprintf('scalp timeseries not consistent with electrode timeseries.\n');
                    error(msg);
                end
            end
            clear TMP;
        end
    end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% surface plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if or(p.elec.plotSurf,p.mesh.plotSurf),
        
        if ~isequal(p.volt.sampleTime,0),
            name = sprintf('Surface: %s @ %8.2f msec',p.volt.file, p.volt.sampleTime);
        else
            name = sprintf('Surface: %s',p.volt.file);
        end
        
        if p.elec.plotSurf,
            
            [p.mesh.current,meshExists] = mesh_check(p,'elec');
            
            if ~meshExists,
                p.mesh.data.meshtype{p.mesh.current}    = 'elec';
                p.mesh.data.vertices{p.mesh.current}    = [Xp,Yp,Zp];
                p.mesh.data.faces{p.mesh.current}       = convhulln([Xp,Yp,Zp]);
                p.mesh.data.lapmat{p.mesh.current}      = [];
                p.mesh.data.lapint{p.mesh.current}      = [];
                p.mesh.data.timeseries{p.mesh.current}  = p.volt.timeArray(:,1);
                p.mesh.data.Cdata{p.mesh.current}       = p.volt.data';
            end
            
            %p = eeg_plot_surf(p);
            %return
            
            H.gui = figure('NumberTitle','off','Name',name,...
                           'PaperType','A4','PaperUnits','centimeters');
            
            set(gca,'Projection','perspective')
            set(gca,'DataAspectRatio',[1 1 1]);
            
            p.elec.patch = patch('vertices',[Xp,Yp,Zp],'faces',convhulln([Xp,Yp,Zp]),...
                            'facevertexCdata',V,'facecolor','interp','EdgeColor','none',...
                            'FaceLighting','phong');
            %patch('vertices',p.elec.data.Vsp,'faces',p.elec.data.Fsp,...
            %      'facevertexCdata',p.elec.data.Psp,'facecolor','interp','EdgeColor','none');
            
            set(gca,'Projection','perspective')
            set(gca,'DataAspectRatio',[1 1 1]);
            axis off tight vis3d
            
            lighting phong
            set(H.gui,'Renderer','zbuffer')
            
            % Create a light above the z axis
            H.light = light('style','infinite','position',[0 0 1000]);
            
            % For a near skin color
            %set(gca,'AmbientLightColor',[.9 .8 .7]);
            
            % MATERIAL([ka kd ks n sc]) sets the ambient/diffuse/specular strength,
            %    specular exponent and specular color reflectance of the objects.
            %p.reflect = material('dull');
            p.reflect{1} = 0.9;
            p.reflect{2} = 0.1;
            p.reflect{3} = 0.0;
            p.reflect{4} = 500;
            p.reflect{5} = 0;
            set(p.elec.patch,'AmbientStrength',p.reflect{1});
            set(p.elec.patch,'DiffuseStrength',p.reflect{2});
            set(p.elec.patch,'SpecularStrength',p.reflect{3});
            set(p.elec.patch,'SpecularExponent',p.reflect{4});
            set(p.elec.patch,'SpecularColorReflectance',p.reflect{5});
            
            if p.elec.plot,
                hold on, plot3(Xp,Yp,Zp,'k.');
            end
            
            set(gca,'Visible','off');
            colormap(p.colorMap.map);
            caxis([p.minimumIntensity p.maximumIntensity]);