[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

    bem_kernel                         C:\BrainStorm_2001\Toolbox\bem_kernel.m
    bst_message_window                 C:\BrainStorm_2001\Toolbox\bst_message_window.m

This function is called by

    bst_headmodeler                    C:\BrainStorm_2001\Toolbox\bst_headmodeler.m

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