[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
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