dipplot.m

含有多种ICA算法的eeglab工具箱

                  'clear editobj tmpobj;' ]);
    h = uicontrol( 'unit', 'normalized', 'position', [0 0.55 .15 .05], 'tag', 'tmp', ...
                  'style', 'pushbutton', 'string', 'Prev', 'callback',  ...
                [ 'editobj = findobj(''parent'', gcf, ''userdata'', ''editor'');' ...
                  'set(editobj, ''string'', num2str(str2num(get(editobj, ''string''))-1));' ...
                  'eval(get(editobj, ''callback''));' ...
                  'clear editobj;' ]);
    h = uicontrol( 'unit', 'normalized', 'position', [0 0.6 .15 .05], 'tag', 'tmp', ...
                  'style', 'pushbutton', 'string', 'Next', 'callback', ...
                [ 'editobj = findobj(''parent'', gcf, ''userdata'', ''editor'');' ...
                  'set(editobj, ''string'', num2str(str2num(get(editobj, ''string''))+1));' ...
                  'dipplot(gcbf);' ...
                  'clear editobj;' ]);
    h = uicontrol( 'unit', 'normalized', 'position', [0 0.65 .15 .05], 'tag', 'tmp', ...
                  'style', 'pushbutton', 'string', 'Keep|Next', 'callback', ...
                [ 'editobj = findobj(''parent'', gcf, ''userdata'', ''editor'');' ...
                  'set(editobj, ''string'', num2str(str2num(get(editobj, ''string''))+1));' ...
                  'tmpobj = get(gcf, ''userdata'');' ...
                  'dipplot(gcbf);' ...
                  'set(tmpobj, ''visible'', ''on'');' ...
                  'clear editobj tmpobj;' ]);
    h = uicontrol( 'unit', 'normalized', 'position', [0 0.7 .15 .05], 'tag', 'tmp', 'userdata', 'rv', ...
                  'style', 'text', 'string', '');
    h = uicontrol( 'unit', 'normalized', 'position', [0 0.75 .15 .05], 'tag', 'tmp', 'userdata', 'editor', ...
                  'style', 'edit', 'string', '1', 'callback', ...
                  [ 'dipplot(gcbf);' ] );
    h = uicontrol( 'unit', 'normalized', 'position', [0 0.8 .15 .05], 'tag', 'tmp', ...
                   'style', 'pushbutton', 'string', 'Plot one', 'callback', ...
             	    [ 'for tmpi = 1:' nbsrc ',' ...
                   '   set(findobj(''parent'', gca, ''tag'', [ ''dipole'' int2str(tmpi) ]), ''visible'', ''off'');' ...
                   'end; clear tmpi;' ...
                   'dipplot(gcbf);' ]);
    h = uicontrol( 'unit', 'normalized', 'position', [0 0.85 .15 .05], 'tag', 'tmp', ...
                   'style', 'pushbutton', 'string', 'Plot All', 'callback', ...
                   [ 'for tmpi = 1:' nbsrc ',' ...
                   '   set(findobj(''parent'', gca, ''tag'', [ ''dipole'' int2str(tmpi) ]), ''visible'', ''on'');' ...
                   'end; clear tmpi;' ]);
    h = uicontrol( 'unit', 'normalized', 'position', [0 0.9 .15 .05], 'tag', 'tmp', ...
                  'style', 'text', 'string', [num2str(length(sources)) ' dipoles:'], 'fontweight', 'bold' );
    h = uicontrol( 'unit', 'normalized', 'position', [0 0.95 .15 .05], 'tag', 'tmp', ...
                  'style', 'text', 'string', '   ' );
    set(gcf, 'userdata', findobj('parent', gca, 'tag', 'dipole1'));

    dat.nbsources  = length(sources);
    set(gca, 'userdata', dat ); % last param=1 for MRI view tight/loose
    set(gcf, 'color', BACKCOLOR);
        
    if strcmp(g.gui, 'off')
        set(findobj('parent', gcf, 'tag', 'tmp'), 'visible', 'off');
    end;
    if strcmp(g.mesh, 'off')
        set(findobj('parent', gca, 'tag', 'mesh'), 'visible', 'off');
    end;
    updatedipplot(gcf);
    
    rotate3d on;
return;

function [xx,yy,zz] = transcoords(x, y, z, TCPARAMS, coreg);
    
    if nargin > 4 & ~isempty(coreg) % custom MRI co-registration
        %                                              rotation    scale     center
        newcoords = transformcoords( [x(:) y(:) z(:)], coreg(5:7), coreg(4), coreg(1:3) );
        xx = newcoords(:,1); xx = reshape(xx, size(x));
        yy = newcoords(:,2); yy = reshape(yy, size(y));
        zz = newcoords(:,3); zz = reshape(zz, size(z));
        return;
    end;
    
    if iscell(TCPARAMS) % coregistration for MRI template
        
        valinion  = TCPARAMS{1};
        valnasion = TCPARAMS{2};
        vallear   = TCPARAMS{3};
        valrear   = TCPARAMS{4};
        valvertex = TCPARAMS{5};
        zoffset   = TCPARAMS{6};
        
        % scale axis
        % ----------
        y = y/2*( valinion(1)  - valnasion(1) );
        x = x/2*( valrear(2)   - vallear(2) );
        z = z*  ( valvertex(3) - vallear(3) - zoffset );
        
        % recenter
        % --------
        yy = y + (vallear (1) + valrear  (1))/2;
        xx = x + (valinion(2) + valnasion(2))/2;
        zz = z + (vallear (3) + valrear  (3))/2 + zoffset;
        return;
    elseif isnumeric(TCPARAMS)
        yy = y*TCPARAMS;
        xx = x*TCPARAMS;
        zz = z*TCPARAMS;        
    end;
    
function sc = plotellipse(sources, ind, nstd, TCPARAMS, coreg);

    for i = 1:length(ind)
        tmpval(1,i) = -sources(ind(i)).posxyz(1);    
        tmpval(2,i) = -sources(ind(i)).posxyz(2);    
        tmpval(3,i) = sources(ind(i)).posxyz(3);
        [tmpval(1,i) tmpval(2,i) tmpval(3,i)] = transcoords(tmpval(1,i), tmpval(2,i), tmpval(3,i), TCPARAMS, coreg);
    end;
    
    % mean and covariance
    C = cov(tmpval');
    M = mean(tmpval,2);
    [U,L] = eig(C);
    
    % For N standard deviations spread of data, the radii of the eliipsoid will
    % be given by N*SQRT(eigenvalues).
    radii = nstd*sqrt(diag(L));
    
    % generate data for "unrotated" ellipsoid
    [xc,yc,zc] = ellipsoid(0,0,0,radii(1),radii(2),radii(3), 10);
    
    % rotate data with orientation matrix U and center M
    a = kron(U(:,1),xc); b = kron(U(:,2),yc); c = kron(U(:,3),zc);
    data = a+b+c;  n = size(data,2);
    x = data(1:n,:)+M(1); y = data(n+1:2*n,:)+M(2); z = data(2*n+1:end,:)+M(3);
    
    % now plot the rotated ellipse
    c = ones(size(z));
    sc = mesh(x,y,z);
    alpha(0.5)
    
function newsrc = convertbesaoldformat(src);
    newsrc = [];
    count = 1;
    countdip = 1;
    if ~isfield(src, 'besaextori'), src(1).besaextori = []; end;
    for index = 1:length(src)
        
        % convert format
        % --------------
        if isempty(src(index).besaextori), src(index).besaextori = 300; end; % 20 mm
        newsrc(count).possph(countdip,:) = [ src(index).besathloc src(index).besaphloc src(index).besaexent];
        newsrc(count).momsph(countdip,:) = [ src(index).besathori src(index).besaphori src(index).besaextori/300];
        
        % copy other fields
        % -----------------
        if isfield(src, 'stdX')
            newsrc(count).stdX = -src(index).stdY;
            newsrc(count).stdY = src(index).stdX;
            newsrc(count).stdZ = src(index).stdZ;
        end;
        if isfield(src, 'rv')
            newsrc(count).rv = src(index).rv;
        end;
        if isfield(src, 'elecrv')
            newsrc(count).rvelec = src(index).elecrv;
        end;
        if isfield(src, 'component')
            newsrc(count).component = src(index).component;
            if index ~= length(src) & src(index).component == src(index+1).component
                countdip = countdip + 1;
            else
                count = count + 1; countdip = 1;
            end;
        else
            count = count + 1; countdip = 1;
        end;
    end; 

function src = computexyzforbesa(src);
    
    for index = 1:length( src )
        for index2 = 1:size( src(index).possph, 1 )

            % compute coordinates
            % -------------------
            postmp = src(index).possph(index2,:);
            momtmp = src(index).momsph(index2,:);
            
            phi      = postmp(1)+90; %% %%%%%%%%%%%%%%% USE BESA COORDINATES %%%%%
            theta    = postmp(2);    %% %%%%%%%%%%%%%%% USE BESA COORDINATES %%%%%
            phiori   = momtmp(1)+90; %% %%%%%%%%%%%% USE BESA COORDINATES %%%%%
            thetaori = momtmp(2);    %% %%%%%%%%%%%% USE BESA COORDINATES %%%%%
            % exentricities are in % of the radius of the head sphere
            [x y z]    = sph2cart(theta/180*pi, phi/180*pi, postmp(3)/100); 
            [xo yo zo] = sph2cart(thetaori/180*pi, phiori/180*pi, momtmp(3)*5); % exentricity scaled for compatibility with DIPFIT
            src(index).posxyz(index2,:) = [-y x z];
            src(index).momxyz(index2,:) = [-yo xo zo];
                    
        end;
     end;
     
% update dipplot (callback call)
% ------------------------------
function updatedipplot(fig)
   
   % find current dipole index and test for authorized range
   % -------------------------------------------------------
   dat     = get(gca, 'userdata');
   editobj = findobj('parent', fig, 'userdata', 'editor');
   tmpnum  = str2num(get(editobj(end), 'string'));
   if tmpnum < 1,             tmpnum = 1;             end;
   if tmpnum > dat.nbsources, tmpnum = dat.nbsources; end;
   set(editobj(end), 'string', num2str(tmpnum));
   
   % hide current dipole, find next dipole and show it
   % -------------------------------------------------
   set(get(gcf, 'userdata'), 'visible', 'off');
   newdip = findobj('parent', gca, 'tag', [ 'dipole' get(editobj(end), 'string')]);
   set(newdip, 'visible', 'on');
   set(gcf, 'userdata', newdip);
   
   % find all dipolar structures
   % ---------------------------
   tmprvobj = findobj('parent', fig, 'userdata', 'rv');
   index   = 1;
   count   = 1;
   for index = 1:length(newdip)
       if isstruct( get(newdip(index), 'userdata') )
           allpos3d(count,:) = getfield(get(newdip(index), 'userdata'), 'pos3d');
           count = count+1;
           foundind = index;
       end;
   end;
 
   % get residual variance
   % ---------------------
   if exist('foundind')
      rv = getfield(get(newdip(foundind), 'userdata'), 'rv');
      set( tmprvobj(end), 'string', rv);
   end
   
   % adapt the MRI to the dipole depth
   % ---------------------------------
   %if ~strcmpi(dat.mode, 'besa') % not besa mode
   delete(findobj('parent', gca, 'tag', 'img'));
   
   tmpdiv = 1;
   if ~dat.axistight
       [xx yy zz] = transcoords(0,0,0, dat.tcparams, dat.coreg);
       indx = minpos(dat.imgcoords{1}-zz);
       indy = minpos(dat.imgcoords{2}-yy);
       indz = minpos(dat.imgcoords{3}-xx);
   else
       if ~dat.cornermri
           indx = minpos(dat.imgcoords{1} - mean(allpos3d(:,3)) + 4/tmpdiv);
           indy = minpos(dat.imgcoords{2} - mean(allpos3d(:,2)) - 4/tmpdiv);
           indz = minpos(dat.imgcoords{3} - mean(allpos3d(:,1)) + 4/tmpdiv);
       else % no need to shift slice if not ploted close to the dipole
           indx = minpos(dat.imgcoords{1} - mean(allpos3d(:,3)));
           indy = minpos(dat.imgcoords{2} - mean(allpos3d(:,2)));
           indz = minpos(dat.imgcoords{3} - mean(allpos3d(:,1)));
       end;
   end;
   plotimgs( dat, [indx indy indz]);
   %end;
   	
% plot images
% -----------
function plotimgs(dat, index);
   
    % loading images
    % --------------
    if strcmpi(dat.mode, 'besa')
        imgt = dat.imgs{1};
        imgc = dat.imgs{2};
        imgs = dat.imgs{3};
    else
        imgt = rot90(squeeze(dat.imgs(:,:,index(1)))); 
        imgc = rot90(squeeze(dat.imgs(:,index(2),:))); 
        imgs = rot90(squeeze(dat.imgs(index(3),:,:))); 
    end;
    if ndims(imgt) == 2, imgt(:,:,3) = imgt; imgt(:,:,2) = imgt(:,:,1); end;
    if ndims(imgc) == 2, imgc(:,:,3) = imgc; imgc(:,:,2) = imgc(:,:,1); end;
    if ndims(imgs) == 2, imgs(:,:,3) = imgs; imgs(:,:,2) = imgs(:,:,1); end;
    
    % computing coordinates for planes
    % --------------------------------    
    wxt = [min(dat.imgcoords{3}) max(dat.imgcoords{3}); min(dat.imgcoords{3}) max(dat.imgcoords{3})];
    wyt = [min(dat.imgcoords{2}) min(dat.imgcoords{2}); max(dat.imgcoords{2}) max(dat.imgcoords{2})];
    wxc = [min(dat.imgcoords{3}) max(dat.imgcoords{3}); min(dat.imgcoords{3}) max(dat.imgcoords{3})];
    wzc = [max(dat.imgcoords{1}) max(dat.imgcoords{1}); min(dat.imgcoords{1}) min(dat.imgcoords{1})];
    wys = [min(dat.imgcoords{2}) max(dat.imgcoords{2}); min(dat.imgcoords{2}) max(dat.imgcoords{2})];
    wzs = [max(dat.imgcoords{1}) max(dat.imgcoords{1}); min(dat.imgcoords{1}) min(dat.imgcoords{1})];
    if dat.axistight & ~dat.cornermri
        wzt = [ 1 1; 1 1]*dat.imgcoords{1}(index(1));
        wyc = [ 1 1; 1 1]*dat.imgcoords{2}(index(2));
        wxs = [ 1 1; 1 1]*dat.imgcoords{3}(index(3));
    else
        wzt = -[ 1 1; 1 1]*dat.maxcoord(1);
        wyc =  [ 1 1; 1 1]*dat.maxcoord(2);
        wxs = -[ 1 1; 1 1]*dat.maxcoord(3);
    end;
    
    % ploting surfaces
    % ----------------
    options = { 'FaceColor','texturemap', 'EdgeColor','none', 'CDataMapping', ...
                'direct','tag','img', 'facelighting', 'none' };
    hold on;
    surface(wxt, wyt, wzt, imgt(end:-1:1,:,:), options{:});
    surface(wxc, wyc, wzc, imgc              , options{:});
    surface(wxs, wys, wzs, imgs              , options{:});
    if strcmpi(dat.drawedges, 'on')
        % removing old edges if any
        delete(findobj( gcf, 'tag', 'edges'));
        if dat.axistight & ~dat.cornermri, col = 'k'; else col = [0.5 0.5 0.5]; end;
        h(1) = line([wxs(1) wxs(2)]', [wys(1) wyc(1)]', [wzt(1) wzt(2)]'); % sagital-transverse
        h(2) = line([wxs(1) wxc(3)]', [wyc(1) wyc(2)]', [wzt(1) wzt(2)]'); % coronal-tranverse
        h(3) = line([wxs(1) wxs(2)]', [wyc(1) wyc(2)]', [wzs(1) wzt(1)]'); % sagital-coronal
        set(h, 'color', col, 'linewidth', 2, 'tag', 'edges');
    end;
        
    %%fill3([-2 -2 2 2], [-2 2 2 -2], wz(:)-1, BACKCOLOR);
    %%fill3([-2 -2 2 2], wy(:)-1, [-2 2 2 -2], BACKCOLOR);
    rotate3d on 

function index = minpos(vals);
	vals(find(vals < 0)) = inf;
	[tmp index] = min(vals);

function scalegca(factor)
    xl = xlim; xf = ( xl(2) - xl(1) ) * factor;
    yl = ylim; yf = ( yl(2) - yl(1) ) * factor;
    zl = zlim; zf = ( zl(2) - zl(1) ) * factor;
    xlim( [ xl(1)-xf xl(2)+xf ]);
    ylim( [ yl(1)-yf yl(2)+yf ]);
    zlim( [ zl(1)-zf zl(2)+zf ]);
    
function color = strcol2real(colorin, colmap)
    if ~iscell(colorin)
        for index = 1:length(colorin)
            color{index} = colmap(colorin(index),:);
        end;
    else
        color = colorin;
        for index = 1:length(colorin)
            if isstr(colorin{index})
                switch colorin{index}
                 case 'r', color{index} = [1 0 0];
                 case 'g', color{index} = [0 1 0];
                 case 'b', color{index} = [0 0 1];
                 case 'c', color{index} = [0 1 1];
                 case 'm', color{index} = [1 0 1];
                 case 'y', color{index} = [1 1 0];
                 case 'k', color{index} = [0 0 0];
                 case 'w', color{index} = [1 1 1];
                 otherwise, error('Unknown color');
                end;
            end;
        end;
    end;