[Master Index]
[Index for Toolbox]
bem_gain
(Toolbox/bem_gain.m in BrainStorm 2.0 (Alpha))
Function Synopsis
bem_gain(bem_grid_mfname,bem_xfer_mfname,ISA,bem_gaingrid_mfname, Verbose)
Help Text
BEM_GAIN - Computes the EEG/MEG forward/gain matrix associated with a set of grid points
function bem_gain(bem_grid_mfname,bem_xfer_mfname,ISA,bem_gaingrid_mfname, Verbose)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
EEG/MEG BEM 3-D GRID FORWARD GAIN CALCULATION UTILITY
This function computes the EEG/MEG forward matrix (kernel) associated with a set of
grid points computed using the utility "gridmaker.m".
INPUTS (Required):
bem_grid_mfname: can be either
- a "*.mat" file name containing pre-computed 3-D grid information
- a structure of specific source locations (3xN array .Loc)and orientations (3xN array .Orient)
bem_xfer_mfname: "*.mat" file containing pre-computed BEM transfer matrix (char string)
ISA: Isolated Skull Approach (0=disable;1=enable) (1 x 1)
bem_gaingrid_mfname: ".mat" output filename used to store final gain matrix and
associated output parameters (char_string)
Verbose: Toggles Verbose mode on (1) / off (0)
INPUTS (Optional):
-None-
OUTPUTS: Complete set of output parameters are stored in user specified
"*.mat" file defined by "bem_gaingrid_mfname.mat"
Other Notes: 3x3 MEG kernel is dotted with sensor orientation to generate more compact
1x3 kernel
External Functions and Files
USC/LANL/CNRS BrainStorm Toolbox
- John J. Ermer 3/7/00
- John J. Ermer 4/28/00 (Added new kernel function to interpolate when using linear basis)
- John J. Ermer 5/05/00 (Added new logic to generate grids for MEG)
- Sylvain Baillet 25/03/02 (Added Verbose mode)
- SB 15-May-2002 Updated Verbose mode with new bst_message_window overwrite features
- SB 07-Mar-2003 Renamed to bem_gain.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cross-Reference Information
This function calls
This function is called by
Listing of function C:\BrainStorm_2001\Toolbox\bem_gain.m
function bem_gain(bem_grid_mfname,bem_xfer_mfname,ISA,bem_gaingrid_mfname, Verbose)
%BEM_GAIN - Computes the EEG/MEG forward/gain matrix associated with a set of grid points
% function bem_gain(bem_grid_mfname,bem_xfer_mfname,ISA,bem_gaingrid_mfname, Verbose)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% EEG/MEG BEM 3-D GRID FORWARD GAIN CALCULATION UTILITY
%
% This function computes the EEG/MEG forward matrix (kernel) associated with a set of
% grid points computed using the utility "gridmaker.m".
%
% INPUTS (Required):
%
% bem_grid_mfname: can be either
% - a "*.mat" file name containing pre-computed 3-D grid information
% - a structure of specific source locations (3xN array .Loc)and orientations (3xN array .Orient)
%
% bem_xfer_mfname: "*.mat" file containing pre-computed BEM transfer matrix (char string)
% ISA: Isolated Skull Approach (0=disable;1=enable) (1 x 1)
% bem_gaingrid_mfname: ".mat" output filename used to store final gain matrix and
% associated output parameters (char_string)
% Verbose: Toggles Verbose mode on (1) / off (0)
% INPUTS (Optional):
% -None-
%
% OUTPUTS: Complete set of output parameters are stored in user specified
% "*.mat" file defined by "bem_gaingrid_mfname.mat"
%
% Other Notes: 3x3 MEG kernel is dotted with sensor orientation to generate more compact
% 1x3 kernel
%
% External Functions and Files
% USC/LANL/CNRS BrainStorm Toolbox
%
% - John J. Ermer 3/7/00
% - John J. Ermer 4/28/00 (Added new kernel function to interpolate when using linear basis)
% - John J. Ermer 5/05/00 (Added new logic to generate grids for MEG)
% - Sylvain Baillet 25/03/02 (Added Verbose mode)
% - SB 15-May-2002 Updated Verbose mode with new bst_message_window overwrite features
% - SB 07-Mar-2003 Renamed to bem_gain.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%<autobegin> ---------------------- 14-Jun-2004 17:09:45 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Forward Modeling
%
% Alphabetical list of external functions (non-Matlab):
% toolbox\bem_kernel.m
% toolbox\bst_message_window.m
%
% At Check-in: $Author: Mosher $ $Revision: 18 $ $Date: 6/14/04 3:37p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 14-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> ------------------------ 14-Jun-2004 17:09:45 -----------------------
%
%%%% THIS PART LOADS AND OPENS INPUT DATA FILES %%%%
%
if ischar(bem_grid_mfname) % Locations of grid points are specified in a file
if ~exist(bem_grid_mfname,'file')
errordlg(sprintf('BEM Grid File %s does not exist on this platform - please re-compute it',bem_grid_mfname))
return
end
load(bem_grid_mfname); % Load Pre-computed3-D Grid Parameters
elseif isstruct(bem_grid_mfname) % Point locations are specified in a structure
Rq_bemgrid = bem_grid_mfname.Loc';
else
errordlg('Location of grid points are passed with wrong argument format',sprintf('Error when calling %s',mfilename))
return
end
if ~exist(bem_xfer_mfname)
errordlg(sprintf('BEM Transfer Matrix File %s does not exist on this platform - please re-compute it',bem_xfer_mfname));
return
end
Xfer = load(bem_xfer_mfname); % Load Pre-computed Transfer Matrices
%
%%%% THIS PART CHECKS IF INPUT PARAMETERS ARE VALID %%%%
%
if (Xfer.modality<1)|(Xfer.modality>2)
error(sprintf('Invalid Modality of %1.0f has been specified!!!',Xfer.modality))
end
if (Xfer.modality==1)& ~isfield(Xfer,'Te')
error(sprintf('EEG Transfer Matrices are not contained in the file: %s',bem_xfer_mfname))
elseif (Xfer.modality==2)& ~isfield(Xfer,'Tm')
error(sprintf('MEG Transfer Matrices are not contained in the file: %s',bem_xfer_mfname))
end
%
%%%% THIS PART LOADS THE APPROPRIATE TRANSFER MATRIX: STANDARD OR ISA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if ISA==0
if Xfer.modality==1
Te_save = Xfer.Te;
Tm_save = []; % Dummy value for kernel routine
else
Tm_save = Xfer.Tm;
Te_save = []; % Dummy value for kernel routine
end
elseif ISA==1
if Xfer.modality==1
Te_save = Xfer.Te_ISA;
Tm_save = [];
else
Tm_save = Xfer.Tm_ISA;
Te_save = [];
end
else
if Verbose
bst_message_window(sprintf('Invalid ISA option chosen!!!'))
end
end
%
clear Xfer.Te Xfer.Tm Xfer.Te_ISA Xfer.Tm_ISA
%
%%%% THIS PART DETERMINES DIMENSION PARAMETERS AND DIVIDES THE INPUT DIPOLE SET INTO A SET OF BLOCKS %%%%
%
if Xfer.modality==1 %EEG
M = size(Xfer.R_eeg,1); % Number of Sensors
S = [];
R_meg = [];
O_meg = [];
elseif Xfer.modality==2 %MEG
M = size(Xfer.R_meg,1);
Xfer.R_eeg = [];
Xfer.P_wts = [];
end
%
P = size(Rq_bemgrid,1); % Number of Dipole Locations
%
BLK_SIZE = 50; % Maximum Number of dipoles to process at a single time
% Set to prevent overstepping memory bounds in Kernel Calculation
blk_tot = floor(P/BLK_SIZE); % Total Number of blocks
blk_rem = rem(P,BLK_SIZE); % Residual Number of Dipoles
%
%%%% THIS PART COMPUTES THE GAIN MATRIX KERNEL (PROCESSING A BLOCK AT A TIME) %%%%
%
t0 = clock; % Start the clock a tickin...
if Verbose
hwait = waitbar(0,'Creating the gain matrix for selected sources. . .');
end
%
if Xfer.modality==1 % EEG
Reeg = Xfer.eegsensor(:,2:4);
Rmeg = [];
GBEM_grid = zeros(M,3*P); % Pre-initialize Output
for nblk = 1:blk_tot
n1 = (nblk-1)*3*BLK_SIZE + 1; % equivalent dipole start index
n2 = nblk*3*BLK_SIZE; % equivalent dipole end index
n3 = (nblk-1)*BLK_SIZE + 1; % dipole start index
n4 = nblk*BLK_SIZE; % dipole end index
[GBEM_grid(:,n1:n2) dummy] = bem_kernel(Rq_bemgrid(n3:n4,:),Xfer.basis,Xfer.test,Xfer.geometry,Xfer.nodes,Xfer.cdv,Te_save,Xfer.P_wts,Tm_save,S);
telap = etime(clock,t0);
trem = telap*(P-n4)/n4;
if Verbose & ~rem(nblk,10)
%bst_message_window('overwrite', sprintf('Completed %.0f of %.0f grid points in %3.2f Sec.',n4,P,telap))
%Expected time to complete: %.2f Sec....',n4,P,telap,trem))
waitbar(nblk/blk_tot, hwait);
end
end
if blk_rem > 0
n1 = blk_tot*3*BLK_SIZE + 1; % equivalent dipole start index (for residual dipoles)
n2 = 3*P; % equivalent dipole end index (for residual dipoles)
n3 = blk_tot*BLK_SIZE + 1; % dipole start index
n4 = P; % dipole end index
[GBEM_grid(:,n1:n2) dummy] = bem_kernel(Rq_bemgrid(n3:P,:),Xfer.basis,Xfer.test,Xfer.geometry,Xfer.nodes,Xfer.cdv,Te_save,Xfer.P_wts,Tm_save,S);
end
elseif Xfer.modality==2 % MEG
Reeg = [];
GBEM_grid = zeros(M,3*P); % Pre-initialize Output
Gtemp = zeros(M,3*BLK_SIZE);
for nblk = 1:blk_tot
n1 = (nblk-1)*3*BLK_SIZE + 1; % equivalent dipole start index
n2 = nblk*3*BLK_SIZE; % equivalent dipole end index
n3 = (nblk-1)*BLK_SIZE + 1; % dipole start index
n4 = nblk*BLK_SIZE; % dipole end index
Ptemp = BLK_SIZE; % number dipoles being processed at current time
[dummy Gtemp] = bem_kernel(Rq_bemgrid(n3:n4,:),Xfer.basis,Xfer.test,Xfer.geometry,Xfer.nodes,Xfer.cdv,Te_save,Xfer.P_wts,Tm_save,Xfer.R_meg);
%
for m=1:M
GBEM_grid(m,n1:n2) = Xfer.O_meg(m,:)*Gtemp(3*m-2:3*m,:); % Apply Sensor Orient to each lead field
end
%
telap = etime(clock,t0);
trem = telap*(P-n4)/n4;
if Verbose & ~rem(nblk,100)
%bst_message_window('overwrite', sprintf('Completed %.0f of %.0f grid points in %3.1f Sec.',n4,P,telap));
%Expected time to complete: %.2f Sec....',n4,P,telap,trem))
waitbar(nblk/blk_tot, hwait);
end
end
if blk_rem > 0
n1 = blk_tot*3*BLK_SIZE + 1; % equivalent dipole start index (for residual dipoles)
n2 = 3*P; % equivalent dipole end index (for residual dipoles)
n3 = blk_tot*BLK_SIZE + 1; % dipole start index
n4 = P; % dipole end index
Ptemp = n4-n3+1; % number dipoles being processed at current time
Gtemp = zeros(M,3*Ptemp);
[dummy, Gtemp] = bem_kernel(Rq_bemgrid(n3:P,:),Xfer.basis,Xfer.test,Xfer.geometry,Xfer.nodes,Xfer.cdv,Te_save,Xfer.P_wts,Tm_save,Xfer.R_meg);
%
for m=1:M
GBEM_grid(m,n1:n2) = Xfer.O_meg(m,:)*Gtemp(3*m-2:3*m,:); % Apply Sensor Orient to each lead field
end
end
end
clear Gtemp
% % If required, apply source orientation
% if isstruct(bem_grid_mfname)
% if ~isempty(bem_grid_mfname.Orient) % Orientations are actually specified
% if Verbose & ~isempty(bem_grid_mfname.Orient)
% bst_message_window('Applying source orientation. . .')
% isrc = 0;
% Gtmp = GBEM_grid;
% clear GBEM_grid
% for src = 1:size(Gtmp,2)/3
% isrc = isrc + 1;
% GBEM_grid(:,src) = Gtmp(:, 3*(isrc-1)+1:3*isrc) * bem_grid_mfname.Orient(:,src);
% end
% end
% end
% end
telap_fwdgaingrid = etime(clock, t0);
if Verbose
close(hwait)
bst_message_window(sprintf('Total time : %3.1f sec.',telap_fwdgaingrid))
end
%
%%%% THIS PART STORES THE FINAL RESULTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% str1 = 'daz az_min az_max del el_min el_max dRc dRf NRc_min Rmin AZ EL Rgrid Rgrid_N';
% str2 = 'Rmax Rbound c0_bemgrid Rq_bemgrid_aer Rq_bemgrid_c0 Rq_bemgrid c0_bemgrid';
% str3 = 'Xfer.R_eeg R_meg O_meg bem_grid_mfname bem_xfer_mfname Xfer.modality GBEM_grid dip_azel_indx';
%
dip_azel_indx = [];
modality = Xfer.modality;
if isfield(Xfer,'R_eeg')
R_eeg = Xfer.R_eeg ;
else
R_eeg = [] ;
end
if isfield(Xfer,'R_meg')
R_meg = Xfer.R_meg ;
O_meg = Xfer.O_meg;
else
R_meg = [];
O_meg = [];
end
str1 = 'R_eeg R_meg O_meg bem_grid_mfname bem_xfer_mfname modality GBEM_grid dip_azel_indx';
%eval(sprintf('save %s %s %s %s',bem_gaingrid_mfname,str1,str2,str3))
eval(sprintf('save %s %s %s %s',bem_gaingrid_mfname,str1))
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