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

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

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