[Master Index] [Index for Toolbox]

# overlapping_sphere

## (Toolbox/overlapping_sphere.m in BrainStorm 2.0 (Alpha))

### Function Synopsis

    Sphere = overlapping_sphere(Channel,SubjectFV, VERBOSE, GRAPHICS)

### Help Text

    OVERLAPPING_SPHERE - Create overlapping spheres for E/MEG channels.
     function Sphere = overlapping_sphere(Channel,SubjectFV, VERBOSE, GRAPHICS)

     Channel is an array of structures.
     Channel.Loc, a matrix of locations of the sensor coil centers, one column per coil
     Channel.Orient, corresponding matrix of orientations, null in the case of EEG
     Channel.Weight, relative weights to each coil. Used here in the least-squares (default is 1).
     Channel.Type, 'MEG','EEG' are the only supported types
     Channel.Comment, unused here.

     SubjectFV.vertices contains the vertices information (3 x Nverts)
     SubjectFV.faces contains the corresponding faces information (Nfaces x 3)
      The faces information is used for determining the surface normals
      and for visualization.

     VERBOSE = 1 (optional), text info.
     GRAPHICS = 1 (optional), graphics displayed

     Sphere is an array of structures
     For each Channel, returns
      Sphere.Center, a 3 x 1 center location and
      Sphere.Radius, the corresponding radius.
      Sphere.Weight is a weighting vector for coloring the vertices to indicate
      the relative surface weights.
      Sphere.Approx is the identified surface weights on the sphere.

      Sphere.Center and Sphere.Radius should be transferred to fields in the
       structure HeadModel.Param

### Cross-Reference Information

This function calls

    bst_message_window                 C:\BrainStorm_2001\Toolbox\bst_message_window.m
    colnorm                            C:\BrainStorm_2001\Toolbox\colnorm.m
    overlapping_sphere_fmins           C:\BrainStorm_2001\Toolbox\overlapping_sphere_fmins.m
    plot3d                             C:\BrainStorm_2001\Toolbox\plot3d.m
    source_grid                        C:\BrainStorm_2001\Toolbox\source_grid.m
    tesselate                          C:\BrainStorm_2001\Toolbox\tesselate.m
    view_surface                       C:\BrainStorm_2001\Toolbox\view_surface.m
    weighting_scalar                   C:\BrainStorm_2001\Toolbox\weighting_scalar.m

This function is called by

    bst_headmodeler                    C:\BrainStorm_2001\Toolbox\bst_headmodeler.m

### Listing of function C:\BrainStorm_2001\Toolbox\overlapping_sphere.m

    function Sphere = overlapping_sphere(Channel,SubjectFV, VERBOSE, GRAPHICS)
    %OVERLAPPING_SPHERE - Create overlapping spheres for E/MEG channels.
    % function Sphere = overlapping_sphere(Channel,SubjectFV, VERBOSE, GRAPHICS)
    %
    % Channel is an array of structures.
    % Channel.Loc, a matrix of locations of the sensor coil centers, one column per coil
    % Channel.Orient, corresponding matrix of orientations, null in the case of EEG
    % Channel.Weight, relative weights to each coil. Used here in the least-squares (default is 1).
    % Channel.Type, 'MEG','EEG' are the only supported types
    % Channel.Comment, unused here.
    %
    % SubjectFV.vertices contains the vertices information (3 x Nverts)
    % SubjectFV.faces contains the corresponding faces information (Nfaces x 3)
    %  The faces information is used for determining the surface normals
    %  and for visualization.
    %
    % VERBOSE = 1 (optional), text info.
    % GRAPHICS = 1 (optional), graphics displayed
    %
    % Sphere is an array of structures
    % For each Channel, returns
    %  Sphere.Center, a 3 x 1 center location and
    %  Sphere.Radius, the corresponding radius.
    %  Sphere.Weight is a weighting vector for coloring the vertices to indicate
    %  the relative surface weights.
    %  Sphere.Approx is the identified surface weights on the sphere.
    %
    %  Sphere.Center and Sphere.Radius should be transferred to fields in the
    %   structure HeadModel.Param

    %<autobegin> ---------------------- 16-Jun-2004 16:50:21 -----------------------
    % --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
    %
    % CATEGORY: Forward Modeling
    %
    % Alphabetical list of external functions (non-Matlab):
    %   toolbox\bst_message_window.m
    %   toolbox\colnorm.m
    %   toolbox\overlapping_sphere_fmins.m
    %   toolbox\plot3d.m
    %   toolbox\source_grid.m
    %   toolbox\tesselate.m
    %   toolbox\view_surface.m
    %   toolbox\weighting_scalar.m
    %
    % At Check-in: $Author: Mosher $  $Revision: 19 $  $Date: 6/16/04 3:16p $
    %
    % This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 16-Jun-2004
    %
    % Principal Investigators and Developers:
    % ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
    %    University of Southern California, Los Angeles, CA
    % ** John C. Mosher, PhD, Biophysics Group,
    %    Los Alamos National Laboratory, Los Alamos, NM
    % ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
    %    CNRS, Hopital de la Salpetriere, Paris, France
    %
    % See BrainStorm website at http://neuroimage.usc.edu for further information.
    %
    % Copyright (c) 2004 BrainStorm by the University of Southern California
    % This software distributed  under the terms of the GNU General Public License
    % as published by the Free Software Foundation. Further details on the GPL
    % license can be found at http://www.gnu.org/copyleft/gpl.html .
    %
    % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
    % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
    % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
    % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
    % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
    %<autoend> ------------------------ 16-Jun-2004 16:50:21 -----------------------

    % /---Script Author--------------------------------------\
    % |                                                      |
    % | *** John C. Mosher, Ph.D.                            |
    % |  Biophysics Group                                    |
    % \------------------------------------------------------/

    % History ----------------------------------------------------
    % Date of creation: Janurary 2001
    % Date of modification: March 2002
    %
    % March 22, 2002 : SB Altered: Script header / response to VERBOSE mode
    %                  / replaced FMINS by FMINSEARCH
    % JCM 08-Sep-2003  Commenting
    % JCM 18-Nov-2003 Replacing .Vertices and .Faces with .vertices and .faces
    %                 in the structure SubjectFV to be
    %                 consistent with the bst_headmodeler call.
    % SB  15-Jun-2004 Minor graphics enhancements
    % ------------------------------------------------------------

    % Calculate the Normals

    if ~exist('VERBOSE','var') | isempty(VERBOSE),
       VERBOSE = 0; % default is silent running
    end

    if ~exist('GRAPHICS','var') | isempty(GRAPHICS),
       GRAPHICS = 0; % default is no display
    end

    %figure(windfind('Junk'))
    % Matlab has the unusual spec of vertices being n x 3, not 3 x n, in patch.
    hp = patch('faces',SubjectFV.faces,'vertices',SubjectFV.vertices','Visible','off');
    Normals = get(hp,'VertexNormals');
    [ignore,Normals] = colnorm(Normals'); % transpose and unit norm of each vertex
    delete(hp)
    %close(gcf)

    Nc = length(Channel); % how many channels to process
    % allocate
    [Sphere(1:Nc)] = deal(struct('Center',[],'Radius',[],'Weight',[],'Approx',[]));

    if VERBOSE
        hwait = waitbar(0,sprintf('Fitting Spheres to %.0f Channels',Nc));
        set(hwait,'units','norm'); % set units to normalized
        hwait_pos = get(hwait,'position');
        set(hwait,'position',[hwait_pos(1) 0.1 hwait_pos(3:4)]); % shift to bottom away from other images
        drawnow
    end

    for i = 1:Nc, % for each channel
        % Calculate the surface weight for the true shape
        if isempty(Channel(i).Weight)
            Channel(i).Weight = 1;
        end

        Sphere(i).Weight = weighting_scalar(Channel(i).Loc,Channel(i).Orient,...
            Channel(i).Weight,SubjectFV.vertices, Normals);

        % initialize the guess of a sphere center as the weighted average of the
        % vertices
        %  tempCenter = mean(SubjectFV.vertices.*[[1;1;1]*Sphere(i).Weight],2); % weighted avg
        tempCenter = mean(SubjectFV.vertices,2);

        % distance from sensor to center is
        tempD = colnorm(Channel(i).Loc(:,1) - tempCenter);

        if(0)
            % initialize radius as average distance to center
            tempR = mean(colnorm(...
                SubjectFV.vertices - tempCenter(:,ones(1,size(SubjectFV.vertices,2)))));
        else
            % 11/5/99 above was nice idea, but sometimes yielded radii larger than sensor distance,
            %  due to neck region. So keep it easy
            tempR = min(.9*tempD,abs(tempD - 0.03)); % come in 3 cm, handle pathological cases
        end

        % grid something comparable
        testL = source_grid(tempD,tempR,1/1000,.3)'; % sparse sampling
        testL = [testL [testL(1:2,:);-testL(3,:)]]; % upper and lower hemispheres

        testL = testL + tempCenter(:,ones(1,size(testL,2)));

        % scan over this grid for weighting scalars, find the best match

        BestGridPt = 1; % assume first
        best_err = overlapping_sphere_fmins([testL(:,1);tempR],...
            Sphere(i).Weight,Channel(i).Loc,Channel(i).Orient,Channel(i).Weight,SubjectFV.vertices);
        for iGrid = 2:size(testL,2), % for each additional grid point
            test_err = overlapping_sphere_fmins([testL(:,iGrid);tempR],...
                Sphere(i).Weight,Channel(i).Loc,Channel(i).Orient,Channel(i).Weight,SubjectFV.vertices);
            if(test_err < best_err),
                BestGridPt = iGrid;
            end
        end

        % set start of search here.
        tempCenter = testL(:,iGrid);
        tempR = mean(colnorm(...
            SubjectFV.vertices - tempCenter(:,ones(1,size(SubjectFV.vertices,2)))));

        if VERBOSE
            bst_message_window({...
                    sprintf('Channel # %d / %d',i,length(Channel)),...
                    sprintf('Initial location is %2.2f %2.2f %2.2f (cm)',100*tempCenter),...
                    sprintf('Initial radius is %2.2f (cm)',100*tempR)})
        end

        X = fminsearch('overlapping_sphere_fmins',[tempCenter;tempR],[],...
            Sphere(i).Weight,Channel(i).Loc,Channel(i).Orient,Channel(i).Weight,SubjectFV.vertices);
        Sphere(i).Center = X(1:3); % the location

        [err,sc,Radius] = overlapping_sphere_fmins(X,...
            Sphere(i).Weight,Channel(i).Loc,Channel(i).Orient,Channel(i).Weight,SubjectFV.vertices);
        Sphere(i).Radius = Radius; % the radius, mapped back depending on the Case Switch.
        Sphere(i).Approx = sc;

        if VERBOSE,

           if GRAPHICS & i == 1,
              fig_os = figure;
              set(fig_os,'Name','Overlapping-Sphere Fitting')
              r = ceil(sqrt(length(Channel))); % Number of rows for following subplot
              c = ceil(sqrt(length(Channel)));
              FigName = get(fig_os,'Name');
           end

           bst_message_window({...
                 sprintf('Final location is %2.2f %2.2f %2.2f (cm)',100*Sphere(i).Center),...
                 sprintf('Final radius is %2.2f (cm)',100*Sphere(i).Radius),' '});

           if GRAPHICS,
              figure(fig_os)
              axx = subplot(r,c,i);
              if i == 1
                  [hf,hs,hl] = view_surface(FigName,...
                      SubjectFV.faces,SubjectFV.vertices,abs(Sphere(i).Weight));
                  set(hs,'facealpha',.8)
                  set(hs,'edgecolor','flat');
                  camlight, lighting gouraud
              else % Spare time by just copying the basic scalp objects
                  C = copyobj([hs,hl],axx);
                  axes(axx)
                  axis off,
                  hs = C(1);
              end
              hold on
              [ignore,ignore,ignore,SVerts,SFaces] = tesselate(Sphere(i).Radius,.15*Sphere(i).Radius,'full');
              SVerts = SVerts'; % columns are vertices
              [ignore,SNormals] = colnorm(SVerts);
              SVerts = SVerts + Sphere(i).Center*ones(1,size(SVerts,2));
              % Channel graph object
              p=patch('vertices',SVerts','faces',SFaces','facecolor',[1 .8 .8],'edgecolor',[1 .5 .5]);
              plot3d(Channel(i).Loc','g+')
              set(p,'Facealpha',.6,'edgecolor','none')

              axis vis3d tight equal
              view(-90,0);
              drawnow

              figure(hwait) % bring to foreground
           end

           waitbar(i/Nc,hwait); % update waitbar
           drawnow

        end

    end

    if(VERBOSE),
        close(hwait);
    end

    return

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Nov 2 16:17:59 2004
Cross-Directory links are: ON