dipplot.m

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

%dipplot normalization
%
%Revision 1.32  2003/06/10 19:04:11  arno
%nothing
%
%Revision 1.31  2003/06/03 16:37:16  arno
%tag images
%
%Revision 1.30  2003/05/30 17:16:22  arno
%nothing
%
%Revision 1.29  2003/05/30 17:09:00  arno
%for index = 1:size(x, 2);
%    dipstruct(index).posxyz = [x(index) y(index) z(index)];
%    %dipstruct(index).posxyz = tmp(index).posxyz;
%    dipstruct(index).momxyz = [1 1 1];
%    dipstruct(index).component = index;
%    dipstruct(index).rv = 0.1;
%end;
%dipplot(dipstruct);
%making xyz output compatible
%
%Revision 1.28  2003/05/30 16:06:27  arno
%nothing
%
%Revision 1.27  2003/05/14 21:36:36  arno
%nothing
%
%Revision 1.26  2003/04/30 18:42:07  arno
%calibrating roughtly the slice selection
%
%Revision 1.25  2003/04/30 16:19:28  arno
%calibrating infants
%
%Revision 1.24  2003/04/30 02:05:24  arno
%changing axis properties for images
%
%Revision 1.23  2003/04/30 01:31:53  arno
%infant option
%
%Revision 1.22  2003/04/23 18:35:11  arno
%allow to plot elipses
%
%Revision 1.21  2003/04/22 21:18:44  arno
%standard dev
%
%Revision 1.20  2003/04/19 01:15:07  arno
%debugging 2 besa dipoles
%
%Revision 1.19  2003/04/19 00:55:53  arno
%correct normalized dipole length
%
%Revision 1.18  2003/04/19 00:46:41  arno
%correcting projection
%
%Revision 1.17  2003/04/19 00:37:43  arno
%changing dipole size for BESA
%
%Revision 1.16  2003/04/11 17:26:45  arno
%accurate plotting in fullMRI
%
%Revision 1.15  2003/04/10 17:37:18  arno
%multi layer MRI plot
%
%Revision 1.14  2003/03/14 17:06:42  arno
%adding error message if plotting non-BESA dipoles on the BESA head model
%
%Revision 1.13  2003/03/14 02:11:33  arno
%automatic scaling for dipfit
%
%Revision 1.12  2003/03/11 23:33:27  arno
%typo
%
%Revision 1.11  2003/03/11 23:27:09  arno
%adding normlen parameter
%
%Revision 1.10  2003/03/11 01:20:32  arno
%updating default besaextori, debuging summary
%
%Revision 1.9  2003/03/07 00:32:53  arno
%debugging textforgui
%
%Revision 1.8  2003/03/06 17:01:10  arno
%textgui -> textforgui
%
%Revision 1.7  2003/03/06 16:50:08  arno
%adding log message
%

function [outsources, XX, YY, ZZ, XO, YO, ZO] = dipplot( sourcesori, varargin )
    
    DEFAULTVIEW = [0 0 1];
        
    if nargin < 1
        help dipplot;
        return;
    end;
        
    % reading and testing arguments
    % -----------------------------
    sources = sourcesori;
    if ~isstruct(sources)
        updatedipplot(sources(1)); 
        % sources countain the figure handler
        return
   end;
    
    %                             key        type       range             default
    g = finputcheck( varargin, { 'color'     ''         []                  [];
                                 'axistight' 'string'   { 'on' 'off' }     'off';
                                 'coreg'     'real'     []                  [];
                                 'drawedges' 'string'   { 'on' 'off' }     'off';
                                 'mesh'      'string'   { 'on' 'off' }     'off';
                                 'gui'       'string'   { 'on' 'off' }     'on';
                                 'summary'   'string'   { 'on' 'off' }     'off';
                                 'view'      'real'     []                 [0 0 1];
                                 'rvrange'   'real'     [0 Inf]             [];
                                 'normlen'   'string'   { 'on' 'off' }     'off';
                                 'num'       'string'   { 'on' 'off' }     'off';
                                 'cornermri' 'string'   { 'on' 'off' }     'off';
                                 'std'       'cell'     []                  {};
                                 'projimg'   'string'   { 'on' 'off' }     'off';
                                 'projcol'   ''         []       [];
                                 'projlines' 'string'   { 'on' 'off' }     'off';
                                 'pointout'  'string'   { 'on' 'off' }     'off';
                                 'dipolesize' 'real'    [0 Inf]             30;
                                 'dipolelength' 'real'  [0 Inf]             1;
                                 'sphere'      'real'   [0 Inf]             1;
                                 'spheres'    'string'  {'on' 'off'}       'off';
                                 'links'       'real'   []                  [];
                                 'image'       { 'string' 'real'} []                 'mri' }, ...
                                                                                    'dipplot');
    if isstr(g), error(g); end;
    g.zoom = 1500;

    % axis image and limits
    % ---------------------
    dat.mode       = g.image;
    dat.maxcoord   = [ 90 100 100 ]; % location of images in 3-D, Z then Y then X
    dat.axistight  = strcmpi(g.axistight, 'on');
    dat.drawedges  = g.drawedges;
    dat.coreg      = g.coreg;
    dat.cornermri  = strcmpi(g.cornermri, 'on');
    radius = 85;
    if isstr(g.image) & strcmpi(g.image, 'besa')
        scaling = 1.05;
        
        % read besa images
        % ----------------
        warning off; imgt = double(imread('besatop.pcx' ))/255; warning on;
        warning off; imgc = double(imread('besarear.pcx'))/255; warning on;
        warning off; imgs = double(imread('besaside.pcx'))/255; warning on;
        dat.imgs        = { imgt imgc imgs };
        
        allcoords1 = ([-1.12 1.12]*size(imgt,2)/200+0.01)*radius; 
        allcoords2 = ([-1.12 1.12]*size(imgt,1)/200+0.08)*radius;
        allcoords3 = [-1.12 1.12]*size(imgs,1)/200*radius; 
        dat.imgcoords = { allcoords3        allcoords2        allcoords1 };

        valinion  = [ 1  0  0 ]*radius;
        valnasion = [-1  0  0 ]*radius;
        vallear   = [ 0 -1  0 ]*radius;
        valrear   = [ 0  1  0 ]*radius;
        valvertex = [ 0  0  1 ]*radius;
        dat.tcparams = { valinion valnasion vallear valrear valvertex 0 };

        %dat.imageaxis   = { -1  1 -1 };
        %dat.imageoffset = { [0.0  0.08 NaN ]  [0.01    NaN -0.01]  [ NaN -0.01   -0.025 ] };
        %dat.imagemult   = { 1.01 1.06 0.96 };
        %dat.axislim     = [-1.2 1.2 -1.2 1.2 -1.2 1.2];
        COLORMESH       = [.5 .5 .5];
        BACKCOLOR       = 'w';
    elseif isstr(g.image)
        fid = fopen('VolumeMNI.bin', 'rb', 'ieee-le');
        if fid == -1
            error('Cannot find MRI data file');
        end;
        V = double(fread(fid, [108 129*129], 'uint8'))/255;
        V = reshape(V, 108, 129, 129);
        fclose(fid);
        
        %V = floatread('VolumeMNI.fdt');
        %load('/home/arno/matlab/MNI/VolumeMNI.mat');
        
        dat.imgs   = V; %smooth3(V,'gaussian', [3 3 3]);
        coordinc   = 2; % 2 mm
        allcoords1 = [0.5:coordinc:size(V,1)*coordinc]-size(V,1)/2*coordinc; 
        allcoords2 = [0.5:coordinc:size(V,2)*coordinc]-size(V,2)/2*coordinc;
        allcoords3 = [0.5:coordinc:size(V,3)*coordinc]-size(V,3)/2*coordinc;
        
        dat.imgcoords = { allcoords3        allcoords2        allcoords1 };
        if strcmpi(g.cornermri, 'on') % make the corner of the MRI volume match
            dat.maxcoord = [max(dat.imgcoords{1}) max(dat.imgcoords{2}) max(dat.imgcoords{3})];
        end;

        COLORMESH = 'w';
        BACKCOLOR = 'k';
        %valinion  = [ 58.5413  -10.5000  -30.8419 ]*2;
        %valnasion = [-56.8767  -10.5000  -30.9566 ]*2;
        %vallear   = [ 0.1040   -59.0000  -30.9000 ]*2;
        %valrear   = [ 0.1040    38.0000  -30.9000 ]*2;
        %valvertex = [ 0.0238   -10.5000   49.8341 ]*2;
        valinion  = [ 52.5413  -10.5000  -30.8419 ]*2;
        valnasion = [-50.8767  -10.5000  -30.9566 ]*2;
        vallear   = [ 0.1040   -51.0000  -30.9000 ]*2;
        valrear   = [ 0.1040    31.0000  -30.9000 ]*2;
        valvertex = [ 0.0238   -10.5000   40.8341 ]*2;
        zoffset   = 27.1190/(27.1190+radius) * (valvertex(3)-vallear(3));
        dat.tcparams = { valinion valnasion vallear valrear valvertex zoffset };
       
        %plotimgs(IMAGESLOC, IMAGESOFFSET, IMAGESMULT, IMAGESAXIS, AXISLIM, [57 85 65]);
        %view(30, 45); axis equal; return;
    else  % custom MRI
        
        V            = -g.image;
        dat.imgs     = -g.image; %smooth3(V,'gaussian', [3 3 3]);
        valinion  = [ 56.5413  0.000  -20.8419 ]*3.5;
        valnasion = [-52.8767  0.000  -20.9566 ]*3.5;
        vallear   = [ -3.1040   -53.0000  -20.9000 ]*3.5;
        valrear   = [ -3.1040    33.0000  -20.9000 ]*3.5;
        valvertex = [ 0.0238   -10.5000   50.8341 ]*3.5;
        zoffset   = 27.1190/(27.1190+radius) * (valvertex(3)-vallear(3));
        dat.tcparams = { valinion valnasion vallear valrear valvertex zoffset };
        dat.coreg    = [];
        coordinc   = 2; % 2 mm
        allcoords1 = [0.5:coordinc:size(V,1)*coordinc]-size(V,1)/2*coordinc; 
        allcoords2 = [0.5:coordinc:size(V,2)*coordinc]-size(V,2)/2*coordinc;
        allcoords3 = [0.5:coordinc:size(V,3)*coordinc]-size(V,3)/2*coordinc;
        dat.imgcoords = { allcoords3        allcoords2        allcoords1 };
        if strcmpi(g.cornermri, 'on') % make the corner of the MRI volume match
            dat.maxcoord = [max(dat.imgcoords{1}) max(dat.imgcoords{2}) max(dat.imgcoords{3})];
        end;

        COLORMESH = 'w';
        BACKCOLOR = 'k';        
    end;

    % point 0
    % -------
    [xx yy zz] = transcoords(0,0,0, dat.tcparams, dat.coreg);
    dat.zeroloc = [ xx yy zz ];
    
    % conversion
    % ----------
    if strcmpi(g.normlen, 'on')
        if isfield(sources, 'besaextori')
            sources = rmfield(sources, 'besaextori'); 
        end;
    end;
    if ~isfield(sources, 'besathloc') & strcmpi(g.image, 'besa') & ~is_sccn
        error(['For copyright reasons, it is not possible to use the BESA ' ...
               'head model to plot non-BESA dipoles']);
    end;
    
    if isfield(sources, 'besathloc')
        sources = convertbesaoldformat(sources);
    end;
    if ~isfield(sources, 'posxyz')
        sources = computexyzforbesa(sources);
    end;        
    if ~isfield(sources, 'component')
        disp('No component indices, making incremental ones...');
        for index = 1:length(sources)
            sources(index).component = index;
        end;
    end;        

    % normalize position to unit sphere
    % ---------------------------------
    maxi = 0;
    for index = 1:length(sources)
        maxi = max(maxi,max(abs(sources(index).posxyz(:))));
    end;
    if maxi > 1.01 & g.sphere == 1
        disp('Non-normalized dipole positions, normalizing by standard head radius 85 mm'); 
        g.sphere = 85;
        fact = 0.1;
    else 
        fact = 1;
    end;
    
    % find non-empty sources
    % ----------------------
    noempt = cellfun('isempty', { sources.posxyz } );
    sources = sources( find(~noempt) );
    
    % transform coordinates
    % ---------------------
    outsources = sources;
    for index = 1:length(sources)
        sources(index).posxyz = sources(index).posxyz/g.sphere;
        tmp = sources(index).posxyz(:,1);
        sources(index).posxyz(:,1) = sources(index).posxyz(:,2);
        sources(index).posxyz(:,2) = -tmp;
        sources(index).momxyz = sources(index).momxyz/g.sphere*0.05*fact;
        tmp = sources(index).momxyz(:,1);
        sources(index).momxyz(:,1) = sources(index).momxyz(:,2);
        sources(index).momxyz(:,2) = -tmp;
        if isfield(sources, 'stdX')
            tmp = sources(index).stdX;
            sources(index).stdX = sources(index).stdY;
            sources(index).stdY = -tmp;
        end;
        if strcmpi(g.normlen, 'on')
            warning off;
            sources(index).momxyz(1,:) = 0.2*sources(index).momxyz(1,:)/ norm(abs(sources(index).momxyz(1,:)));
            if size(sources(index).momxyz,1) > 1 & sources(index).momxyz(1) ~= 0
                sources(index).momxyz(2,:) = 0.2*sources(index).momxyz(2,:)/ norm(abs(sources(index).momxyz(2,:)));
            end;
            warning on;
        end;
    end;
    
    % remove sources with out of bound Residual variance
    % --------------------------------------------------
    if isfield(sources, 'rv') & ~isempty(g.rvrange)
        for index = length(sources):-1:1
            if sources(index).rv < g.rvrange(1)/100 | sources(index).rv > g.rvrange(2)/100
                sources(index) = [];
            end;
        end;
    end;
    
    % color array
    % -----------
    if isempty(g.color)
        g.color = { 'g' 'b' 'r' 'm' 'c' 'y' };
        if strcmp(BACKCOLOR, 'w'), g.color = { g.color{:} 'k' }; end;
    end;
    g.color = g.color(mod(0:length(sources)-1, length(g.color)) +1);
    if ~isempty(g.color)
        g.color = strcol2real( g.color, jet(64) );
    end;
    if ~isempty(g.projcol)
        g.projcol = strcol2real( g.projcol, jet(64) );
        g.projcol = g.projcol(mod(0:length(sources)-1, length(g.projcol)) +1);        
    else
        g.projcol = g.color;
        for index = 1:length(g.color)
            g.projcol{index} =  g.projcol{index}/2;
        end;
    end;
    
    % build summarized figure
    % -----------------------
    if strcmp(g.summary, 'on')
        figure;
        options = { 'gui', 'off', 'dipolesize', g.dipolesize/1.5,'dipolelength', g.dipolelength, 'sphere', g.sphere ...
                    'color', g.color, 'mesh', g.mesh, 'num', g.num, 'image', g.image 'normlen' g.normlen };
        axes('position', [0 0 0.5 0.5]);  dipplot(sourcesori, 'view', [1 0 0] , options{:}); axis off; if strcmpi(g.image, 'besa'), scalegca(0.1); end;