[Master Index] [Index for Toolbox]

# source_grids

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

### Function Synopsis

    HeadModel = source_grids(HeadModel,Channel,SubjectFV,GUI);

### Help Text

    SOURCE_GRIDS - Generate source grids for parametric model searches
    function HeadModel = source_grids(HeadModel,Channel,SubjectFV,GUI);

     HeadModel and Channel are the standard structures:
     HeadModel:
      .Function
      .Param, struct, e.g. 'center','radius'
      .GridLoc, RETURNED: cell array generated by this program
      .SearchTess, unused, maybe dropped
      .Gain, RETURNED: cell array one-to-one with GridLoc
      .GridName
      .GridLoc
      .ImageGridOrient
      .Gain
      .GainCovar

     Because of the size of the Image Grids, I suggest that .Image* fields NOT
      be sent along in the calling structure, but added in later.

     Channel:
      .Loc
      .Orient
      .Weight
      .Type
      .Comment

     SubjectFV is the faces and vertices from the SubjectTess structure. Note the
      structures specific to Matlab, not BrainStorm: lowercase names and row orientation.
     SubjectFV.faces is n x 3, used only to reducepatch the vertices to around 250.
     SubjectFV.vertices is 3 x m, used to bound the search grid, should
      be a smooth inner-skull-like surface inside of which the points will be generated.

     GUI : a structure passing several parameters to the source_grids
     .VALIDORDER = Source orders Handled by the forward models.
           As of 03/13/00 : [-1 0 1] in MEG and [-1] for EEG
           Initially non-synchronous sources only
     .SPACING = ; initial grid spacing in mm for the SearchGrid
     .VERBOSE : 1 (res. 0) turns on (res. off) the verbose mode
     .MEG and .EEG = two binary flags indicating whether the grid gain needs to be computed (1) for either modality or not (0)

### Cross-Reference Information

This function calls

    bst_message_window                 C:\BrainStorm_2001\Toolbox\bst_message_window.m
    closest                            C:\BrainStorm_2001\Toolbox\closest.m
    colnorm                            C:\BrainStorm_2001\Toolbox\colnorm.m
    good_channel                       C:\BrainStorm_2001\Toolbox\good_channel.m
    meg_bem                            C:\BrainStorm_2001\Toolbox\meg_bem.m
    rownorm                            C:\BrainStorm_2001\Toolbox\rownorm.m
    subcorr_gui                        C:\BrainStorm_2001\Toolbox\subcorr_gui.m

This function is called by

    bst_headmodeler                    C:\BrainStorm_2001\Toolbox\bst_headmodeler.m

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

    function HeadModel = source_grids(HeadModel,Channel,SubjectFV,GUI);
    %SOURCE_GRIDS - Generate source grids for parametric model searches
    %function HeadModel = source_grids(HeadModel,Channel,SubjectFV,GUI);
    %
    % HeadModel and Channel are the standard structures:
    % HeadModel:
    %  .Function
    %  .Param, struct, e.g. 'center','radius'
    %  .GridLoc, RETURNED: cell array generated by this program
    %  .SearchTess, unused, maybe dropped
    %  .Gain, RETURNED: cell array one-to-one with GridLoc
    %  .GridName
    %  .GridLoc
    %  .ImageGridOrient
    %  .Gain
    %  .GainCovar
    %
    % Because of the size of the Image Grids, I suggest that .Image* fields NOT
    %  be sent along in the calling structure, but added in later.
    %
    % Channel:
    %  .Loc
    %  .Orient
    %  .Weight
    %  .Type
    %  .Comment
    %
    % SubjectFV is the faces and vertices from the SubjectTess structure. Note the
    %  structures specific to Matlab, not BrainStorm: lowercase names and row orientation.
    % SubjectFV.faces is n x 3, used only to reducepatch the vertices to around 250.
    % SubjectFV.vertices is 3 x m, used to bound the search grid, should
    %  be a smooth inner-skull-like surface inside of which the points will be generated.
    %
    % GUI : a structure passing several parameters to the source_grids
    % .VALIDORDER = Source orders Handled by the forward models.
    %       As of 03/13/00 : [-1 0 1] in MEG and [-1] for EEG
    %       Initially non-synchronous sources only
    % .SPACING = ; initial grid spacing in mm for the SearchGrid
    % .VERBOSE : 1 (res. 0) turns on (res. off) the verbose mode
    % .MEG and .EEG = two binary flags indicating whether the grid gain needs to be computed (1) for either modality or not (0)
    %

    %<autobegin> ---------------------- 26-May-2004 11:34:31 -----------------------
    % --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
    %
    % CATEGORY: Forward Modeling
    %
    % Alphabetical list of external functions (non-Matlab):
    %   toolbox\bst_message_window.m
    %   toolbox\closest.m
    %   toolbox\colnorm.m
    %   toolbox\dip.m
    %   toolbox\good_channel.m
    %   toolbox\rownorm.m
    %   toolbox\subcorr_gui.m
    %
    % At Check-in: $Author: Mosher $  $Revision: 19 $  $Date: 5/26/04 10:02a $
    %
    % This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 24-May-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> ------------------------ 26-May-2004 11:34:31 -----------------------

    % /---Script Authors--------------------------------------\
    % |                                                      |
    % | *** John C. Mosher, Ph.D.                            |
    % |  Biophysics Group                                    |
    % |                                                      |
    % | *** Sylvain Baillet Ph.D.                            |
    % | Cognitive Neuroscience & Brain Imaging Laboratory    |
    % | CNRS UPR640 - LENA                                   |
    % | Hopital de la Salpetriere, Paris, France             |
    % | sylvain.baillet@chups.jussieu.fr                     |
    % |                                                      |
    % \------------------------------------------------------/

    % Script History -----------------------------------------------------------------------------------------------------------
    % Date of creation: October, 25 1999
    %
    % Nov 2001: Checked sequential computation of EEG and MEG head models when requested
    % SB March 2002 :  - Added VERBOSE mode
    %                  - Can handle simulateous creation of MEG and EEG GridLoc and SearchGridGain
    % SB 24-Apr-2002 : - Script header platic surgery
    %                  - Moved back to .Faces and .Vertices field names for SubjectFV
    % SB 26-Apr-2002 : - Moved to BsT MMII new HeadModel filednames (e.g. .GridName)
    %                  - Removed computation of Gains for CD or CME pairs : will sample them from original CD and CME Gains (like we already did actually)
    %                  - Added the GUI.MEG and GUI.EEG flags
    % SB 01-Jul-2002 : - Added error detection when EEG channels are used for EOG or ECoG.
    % SB 02-Jul-2002 : - Fixed bug when gridding for CME only
    % SB 12-Sep-2002 : - SubjectFV is a tessellation structure that depicts the inner-skull volume
    %                    where the search volumic grid must be contained
    %                    Former version used the scalp surface which triggered errors in the 3-shell Legendre gain matrix computation
    %                    because it checked for sources outside the inner sphere.
    % SB 19-Nov-2002 : - Minor display upgrade of display
    %                    Updated management of EEG reference.
    % JCM 20-Nov-2002    Minor fix of duplicate autocomments blocks, some comment formatting
    % JCM 18-Nov-20003   Fixing problem of vertice being m x 3 or 3 x m. Force to be 3 x m.
    % ------------------------------------------------------------------

    SPACING = GUI.SPACING/1000; % in meters
    VALIDORDER = GUI.VALIDORDER;
    SourceOrderString = {'Current Dipole',...
          'Unvalid Order (was Magnetic Dipole Moment)',...
          'Current Multipole Expansion',...
          'Current Dipole Pairs',...
          'Unvalid Order (was Magnetic Dipole Moment PAIRS)',...
          'Current Multipole Expansion Pairs'};

    % Local hardwired options
    SENSOR_PROXIMITY = 25/1000; % grid can't be closer than this to a sensor coil

    % pick a pruning technique
    %PRUNE = 'Power'; % based on the remnant power in the projected gain matrices
    PRUNE = 'Subspace'; % based on the correlation of adjacent points.
    VERBOSE = 0;
    % Pruning Options
    LOCAL = SPACING*2; % local sphere to check pruning options
    CORR = 0.95; % grid points must be less correlated than this to remain good

    [m,n] = size(SubjectFV.vertices);

    if m ~= 3,
       SubjectFV.vertices = SubjectFV.vertices'; % should be 3 x n, each column a vector
    end

    % MEG & EEG respective channel indices
    % Find MEG and EEG CHANndx
    ChannelFlag = ones(size(Channel));
    MEGndx = good_channel(Channel,ChannelFlag,'MEG');
    EEGndx = good_channel(Channel,ChannelFlag,'EEG');

    HeadModel.GridLoc = cell(1,length(VALIDORDER));
    HeadModel.Gain = cell(1,length(VALIDORDER));

    % Check valid source orders depeding on the requests
    if ~GUI.MEG | (GUI.MEG & strcmp(HeadModel.Function{MEGndx(1)},'meg_bem')) % 'doing MEG BEM or EEG only -> Authorize current dipole models
       inotvalid = find(...
          VALIDORDER == 0 | ...
          VALIDORDER == 1 | ...
          VALIDORDER == 3 | ...
          VALIDORDER == 4);

       VALIDORDER(inotvalid) = []; % remove all unvalid source configurations

    end

    for ORDER = VALIDORDER, %the source model orders
       iOrder = ORDER + 2; % associated indexer, starting with 1.

       if 0,%GUI.VERBOSE
          bst_message_window({sprintf('Creating %s gain matrices for volume source scan ', SourceOrderString{iOrder})}),
       end

       switch ORDER % now translate order to model orders
       case {-1,0,1}
          % order remains the same
          SYNC = 1; % single source
       case {2,3,4}
          ORDER = ORDER - 3; % revert to true order
          SYNC = 2; % dual synchronous source
       end

       % We anticipate that single dipoles need the densest grid. Each subsequent
       %  order uses the previous grid.

       if ORDER==VALIDORDER(1), % first time through
          % first generate a plaid grid bounded by the subject vertices

          % only need about 250 vertices to box in the grid.
          if size(SubjectFV.faces,1)>500
             tempFV = SubjectFV;
             tempFV.vertices = tempFV.vertices'; % so that they are n x 3
             NFV = reducepatch(tempFV,500); % number of faces
             GridBox = NFV.vertices'; % 3 x ~500
          else
             GridBox = SubjectFV.vertices;
          end

          mx = max(GridBox,[],2);
          mn = min(GridBox,[],2);

          [X,Y,Z] = meshgrid([mn(1):SPACING:mx(1)],...
             [mn(2):SPACING:mx(2)],[mn(3):SPACING:mx(3)]);

          L = [X(:) Y(:) Z(:)]';
          tempLn = rownorm(L);
          L(find(~tempLn),:) = []; % delete points at the origin as degenerate in MEG

          % prune sources outside the the GridBox
          % We find the vertex closest to the grid point. If the radius of that
          %  grid point is greater than the vertex, we mark it bad.
          % Inelegant 3-d measure of "inside", could use better.
          % "radius" is defined to a central mean of the vertices.
          mnVert = mean(SubjectFV.vertices,2); % mean of the original vertices.

          Good = ones(1,size(L,2));
          if GUI.VERBOSE
             hwait = waitbar(0,sprintf('Sifting %.0f dipoles plaid at %.1f mm',...
                size(L,2),SPACING*1000));
          end

          for i = 1:size(L,2),
             if(~rem(i,100)),
                if GUI.VERBOSE, waitbar(i/size(L,2),hwait); end
             end
             pt2 = closest(L(:,i)',GridBox'); % old routine using rows of data not cols.
             if(norm(L(:,i)-mnVert) > 0.99*norm(pt2'-mnVert)),
                Good(i) = 0; % mark it bad
             end
          end
          if GUI.VERBOSE
             close(hwait), end

          if GUI.MEG & GUI.EEG
             ChanNdx = MEGndx;
          elseif GUI.MEG
             ChanNdx = MEGndx;
          elseif GUI.EEG
             ChanNdx = EEGndx;
          end

          % should also prune out stuff adjacent sensors.
          if GUI.VERBOSE
             hwait = waitbar(0,sprintf('Checking %.0f Sensors for %.1f mm Proximity',...
                length(Channel(ChanNdx)),SENSOR_PROXIMITY*1000));
          end

          for i =1:length(Channel(ChanNdx)),
             if(~rem(i,50)),
                if GUI.VERBOSE, waitbar(i/length(Channel(ChanNdx)),hwait); end
             end
             for j = 1:size(Channel(ChanNdx(1)).Loc,2), % each coil
                DIST = colnorm(Channel(ChanNdx(i)).Loc(:,j+zeros(1,size(L,2))) - L);
                ndx= find(DIST < SENSOR_PROXIMITY); % don't get closer than this
                Good(ndx) = 0; % mark bad
             end
          end
          if GUI.VERBOSE, close(hwait), end
          Gndx = find(Good); % the good grid points
          L = L(:,Gndx); % downselect
       else
          % next order, use the previous grid to start
          L = HeadModel.GridLoc{iOrder-2}; % previous grid (Synchronous sources)
       end

       G = cell(1,2);

       if xor(GUI.MEG,GUI.EEG)

          if GUI.MEG
             ChanNdx = MEGndx;
             sstring = 'MEG';
          else
             ChanNdx = EEGndx;
             sstring = 'EEG';
          end

          if GUI.VERBOSE
             bst_message_window(sprintf('Creating %s gain matrix for %.0f sensors and %.0f sources...',...
                sstring, length(ChanNdx),size(L,2)));
          end

          %G = feval(HeadModel.Function{ChanNdx(1)},L,Channel(ChanNdx),HeadModel.Param(ChanNdx),ORDER); % same for synchronous sources as well
          G = feval(HeadModel.Function{ChanNdx(1)},L,Channel,HeadModel.Param,ORDER); % same for synchronous sources as well

       else % Both MEG and EEG are requested
          if GUI.VERBOSE
             bst_message_window(sprintf('Creating MEG gain matrix first, for %.0f sensors and %.0f sources...',...
                length(MEGndx),size(L,2)));
          end

          %G = feval(HeadModel.Function{MEGndx(1)},L,Channel(MEGndx),HeadModel.Param(MEGndx),ORDER); % same for synchronous sources as well
          G = feval(HeadModel.Function{MEGndx(1)},L,Channel,HeadModel.Param,ORDER); % same for synchronous sources as well

       end

       DIM = size(G,2)/size(L,2);% dimension of each source

       % done generating initial grid and gain matrix

       % now prune the grid of redundant source points

       Good = ones(1,size(L,2)); % every grid point initialized as good
       Goodndx = reshape([1:size(G,2)]',DIM,size(G,2)/DIM); % indexer into gain matrix, by columns
       % therefore Goodndx(DIM,:)/DIM is an indexer into Good.

       PpG = zeros(size(G)); % zero out all locations
       if(strcmp(PRUNE,'Power')),
          % need a knee calculation, where NOS is the std dev of the white noise
          KNEE = NOS*sqrt(size(G,1)); % sum of the total sensor noise
       end
       if GUI.VERBOSE
          hwait = waitbar(0,sprintf('Processing %.0f sources using %s',...
             length(Good),PRUNE)); end

       for i = 1:length(Good),
          NumGoodSources = sum(Good); % how many good sources are there still left
          if(Good(i)), % if it's still good, i.e. hasn't been deleted

             dndx = [-(DIM-1):0]+i*DIM; % indexer to gain matrix

             % Check only those sources in geometric proximity to the source.
             %  Done for speed considerations.
             DIST = colnorm(L(:,i+zeros(1,NumGoodSources)) - L(:,find(Good)));
             localndx = find(DIST <= LOCAL);
             if 0%GUI.VERBOSE
                bst_message_window(sprintf('Processing %.0f with %.0f local out of %.0f good sources. . .',...
                   i,length(localndx),NumGoodSources));            end

             switch PRUNE
             case 'Power'

                Ug = orth(G(:,dndx)); % orthogonalize the source
                % project source matrices away from the present point.
                % but only bother with the good points. The other columns will be left alone.
                %PpG = G - Ug*(Ug'*G);
                PpG(:,Goodndx(:,localndx)) = G(:,Goodndx(:,localndx)) - ...
                   Ug*(Ug'*G(:,Goodndx(:,localndx)));
             case 'Subspace'
                % do nothing
             end

             if GUI.VERBOSE & (~rem(i,100))
                waitbar(i/length(Good),hwait); end

             for j = localndx(end:-1:1), % for all remaining good sources, in reverse search order
                Source_Num = Goodndx(DIM,j)/DIM; % the source number
                switch PRUNE
                case 'Power'
                   Strongest_Source = svd(PpG(:,Goodndx(:,j))); % decompose the projected matrix
                   Strongest_Source = Strongest_Source(1); % the strongest
                   Strongest_Source = Strongest_Source * DIP; % scaled by anticipated dipole strength

                   % now standard deviation of the noise is NOS, the "knee" in the equation.
                   % if the Strongest_Source is below the knee, delete it.

                   if(Strongest_Source < KNEE), % weaker than noise
                      if(Source_Num ~= i), % but don't delete the source itself
                         if 0%GUI.VERBOSE
                            bst_message_window(sprintf('Deleting Source %.0f',Source_Num));
                         end
                         Good(Source_Num) = 0; % mark it deleted
                         Goodndx(:,j) = []; % remove this column
                         % because of reverse sequencing, has no effect on the lower indices in the for loop
                      end
                   end
                case 'Subspace'
                   nextCorr= subcorr_gui(G(:,dndx),G(:,Goodndx(:,j)),...
                      struct('REG','Tikhonov','Tikhonov',100));
                   if isempty(nextCorr) % Occurs when G is full of zeros. Happens for electrodes wrongly involved in forward modeling though they were measuring non-EEG data (like EOG or ECoG)
                      errordlg('Please make sure sensor coordinates are all correct and non-zero','Forward modeling aborted')
                      return
                   end

                   if(nextCorr(end) > CORR), % all corrs too close in correlations
                      if(Source_Num ~= i), % but don't delete the source itself
                         if 0%GUI.VERBOSE,
                            bst_message_window(sprintf('Deleting Source %.0f',Source_Num));
                         end
                         Good(Source_Num) = 0; % mark it deleted
                         Goodndx(:,j) = []; % remove this column
                         % because of reverse sequencing, has no effect on the lower indices in the for loop
                      end
                   end
                end
             end % for loop, testing adequacy of sources
          end % if this source is still good

       end % for all source combinations

       if(GUI.VERBOSE),
          close(hwait), end

       if GUI.EEG
          % Detect EEG Reference
          RefNdx = good_channel(Channel,[],'EEG REF');
          if isempty(RefNdx) % Average reference
             RefNdx = 0;
          end
       end

       if GUI.MEG
          ChanNdx = MEGndx;
       elseif GUI.EEG
          ChanNdx = EEGndx;
       end

       HeadModel.Gain{iOrder} = NaN*ones(length(Channel),length(Goodndx(:)));

       if GUI.VERBOSE
          bst_message_window(sprintf('Decimating source grid to %.0f most significant sources',...
             length(find(Good)))); end

       HeadModel.Gain{iOrder}(ChanNdx,:) = G(:,Goodndx(:)); clear G Goodndx
       HeadModel.GridLoc{iOrder} = L(:,find(Good)); clear L Good

       % Now compute EEG gain matrix from MEG's source grid locations if MEG and EEG are available
       if GUI.MEG & GUI.EEG
          if GUI.VERBOSE
             bst_message_window(sprintf('Now creating EEG gain matrix for %.0f sensors and %.0f sources',...
                length(EEGndx),size(HeadModel.GridLoc{iOrder} ,2))); end

          if iOrder == 1 | iOrder == 4 % Current dipole source models
             G = feval(HeadModel.Function{EEGndx(1)},HeadModel.GridLoc{iOrder},Channel,HeadModel.Param,ORDER); % same for synchronous sources as well
          else
             G = NaN*ones(size(HeadModel.Gain{iOrder}(EEGndx,:)));
          end

          HeadModel.Gain{iOrder}(EEGndx,:) = G;%G(EEGndx,:);

       end

       clear G

       if 0, %GUI.VERBOSE
          bst_message_window({sprintf('Creating %s gain matrices for volume source scan -> DONE', SourceOrderString{iOrder}),' '})
       end

    end % per source model order

    % % for now, let's handle the synchronous sources somewhat crudely. I'm gonna
    % %  just use the multipolar grid, and let rapmusic_gui sample pairs from that
    %
    % for i = VALIDORDER
    %     iOrder = i+2;
    %     HeadModel.GridLoc{iOrder+3} = HeadModel.GridLoc{iOrder};
    %     HeadModel.Gain{iOrder+3} = HeadModel.Gain{iOrder};
    % end

    return

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Oct 12 12:05:14 2004
Cross-Directory links are: ON