[Master Index] [Index for Toolbox]

# bem_xfer

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

### Function Synopsis

    [varargout] = bem_xfer(R_eeg,R_meg,O_meg,vertices,faces,sigma,mode, ...

### Help Text

    BEM_XFER - Computation of the MEG and/or EEG BEM transfer matrices
     function [varargout] = bem_xfer(R_eeg,R_meg,O_meg,vertices,faces,sigma,mode, ...
         basis_opt,test_opt,ISA,fn_eeg,fn_meg,NVertMax,Verbose)
     function [Te_ISA,Te,Tm_ISA,Tm,nfv] = bem_xfer(R_eeg,R_meg,O_meg,vertices,faces,sigma,mode, ...
         basis_opt,test_opt,ISA,fn_eeg,fn_meg,NVertMax,Verbose)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     EEG/MEG FORWARD MODEL USING BOUNDARY ELEMENT METHODS
                          - TRANSFER MATRIX CALCULATION (bem_xfer.m)

     This function computes the "transfer matrix" associated with the voltage
     potential  forward gain matrix for an array of EEG electrodes located on the
     outermost layer of a single/multilayer surface,
                                  -AND/OR-
     the radial magnetic field forward gain matrix for an array of MEG sensors
     located on the  outermost layer of a single/multilayer surface.

     Surface data for each layer is defined via a user specified tesselated grid.
     Each region  of the multilayer surface is assumed to be concentric with
     isotropic conductivity.   EEG sensors are assumed to be located on the surface
     of the outermost layer at one of the surface tesselation points. Should
     specified EEG sensor points not coincide with a point   on the tesselated
     outer surface layer grid, the sensor location will be quantized to the
     closest outer surface tesselation point. Dipole generator(s) are assumed to be
     interior  to the innermost "core" layer

     INPUTS (Required)
           R_eeg  : EEG sensor locations on scalp (meters)         (Meeg x 3)
                      (insert dummy value if mode = 2)
           R_meg  : MEG sensor locations on scalp (meters)         (Mmeg x 3)
                      (insert dummy value if mode = 1)
           O_meg  : MEG sensor orientations (unit-vector)          (Mmeg x 3)
                      (insert dummy value if mode = 1)
        vertices  : Tessalated Surface layers from INNERMOST TO OUTERMOST.
                    Packed as Cell Array where index 1 is innermost surface  (1 x NL)-Cell Array
                    Individual entries are position row vectors in PCS
           faces  : Tesselated Surface Node connections for each triangle    (1 x NL)-Cell Array
                    Packed as Cell Array where index 1 is innermost surface
                    Individual entries are indexes to "vertices"
           sigma  : conductivity from INNERMOST to OUTERMOST        (NL x 1)
           mode   : 1,  compute EEG only
                    2,  compute MEG only
                    3,  compute both EEG and MEG
       basis_opt  : 0, constant basis
                    1, linear
       test_opt   : 0, collocation,
                    1, Galerkin
         ISA      : 0,  Inhibit Isolated Skull Approach
                  : 1,  Enable Isolated Skull Approach

        fn_eeg    : EEG transfer matrix filename                   Char String
                       - EEG Transfer matrices are computed and then saved
                         in this "*.mat" file

        fn_meg    : MEG transfer matrix filename                   Char String
                       - MEG Transfer matrices are computed and then saved
                         in this "*.mat" file

             WHERE: M = # of sensors; P = # of dipoles; NL = # of sphere layers
                   Nn = # of Surface Tesselation Nodes; Nt = # of Surf Tess Triangles

     INPUTS (Optional):
        NVertMax : Maximum Number of vertices for Surface Tesselation  (for a single layer)
                       - Surface is re-tesselated if # points exceeds max
                       - Default Maximum Number of Triangle Faces is 2292    (1 x NL)

        Verbose : Toggles Verbose mode on/off (on = 1);

     OUTPUTS:
        nfv    : a surface tessellation structure with fields .vertices and
                 .faces of all envelopes used in BEM computation, once
                 downsampling, separation and embedding check have been
                 completed.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

### Cross-Reference Information

This function calls

    bst_message_window                 C:\BrainStorm_2001\Toolbox\bst_message_window.m
    crossprod                          C:\BrainStorm_2001\Toolbox\crossprod.m
    dotprod                            C:\BrainStorm_2001\Toolbox\dotprod.m
    get_user_directory                 C:\BrainStorm_2001\Toolbox\get_user_directory.m
    hrow_linear                        C:\BrainStorm_2001\Toolbox\hrow_linear.m
    inpolyhd                           C:\BrainStorm_2001\PublicToolbox\OtherTools\inpolyhd.m
    kl_galerkin                        C:\BrainStorm_2001\Toolbox\kl_galerkin.m
    kl_gram                            C:\BrainStorm_2001\Toolbox\kl_gram.m
    makeuswait                         C:\BrainStorm_2001\Toolbox\makeuswait.m
    ray_intersect                      C:\BrainStorm_2001\Toolbox\ray_intersect.m
    rownorm                            C:\BrainStorm_2001\Toolbox\rownorm.m
    separation                         C:\BrainStorm_2001\Toolbox\separation.m
    solid_angle2                       C:\BrainStorm_2001\Toolbox\solid_angle2.m
    tessellation_outwards              C:\BrainStorm_2001\Toolbox\tessellation_outwards.m
    trans_matrix                       C:\BrainStorm_2001\Toolbox\trans_matrix.m
    view_surface                       C:\BrainStorm_2001\Toolbox\view_surface.m

This function is called by

    bst_headmodeler                    C:\BrainStorm_2001\Toolbox\bst_headmodeler.m

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

    function [varargout] = bem_xfer(R_eeg,R_meg,O_meg,vertices,faces,sigma,mode, ...
        basis_opt,test_opt,ISA,fn_eeg,fn_meg,NVertMax,Verbose)
    %BEM_XFER - Computation of the MEG and/or EEG BEM transfer matrices
    % function [varargout] = bem_xfer(R_eeg,R_meg,O_meg,vertices,faces,sigma,mode, ...
    %     basis_opt,test_opt,ISA,fn_eeg,fn_meg,NVertMax,Verbose)
    % function [Te_ISA,Te,Tm_ISA,Tm,nfv] = bem_xfer(R_eeg,R_meg,O_meg,vertices,faces,sigma,mode, ...
    %     basis_opt,test_opt,ISA,fn_eeg,fn_meg,NVertMax,Verbose)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % EEG/MEG FORWARD MODEL USING BOUNDARY ELEMENT METHODS
    %                      - TRANSFER MATRIX CALCULATION (bem_xfer.m)
    %
    % This function computes the "transfer matrix" associated with the voltage
    % potential  forward gain matrix for an array of EEG electrodes located on the
    % outermost layer of a single/multilayer surface,
    %                              -AND/OR-
    % the radial magnetic field forward gain matrix for an array of MEG sensors
    % located on the  outermost layer of a single/multilayer surface.
    %
    % Surface data for each layer is defined via a user specified tesselated grid.
    % Each region  of the multilayer surface is assumed to be concentric with
    % isotropic conductivity.   EEG sensors are assumed to be located on the surface
    % of the outermost layer at one of the surface tesselation points. Should
    % specified EEG sensor points not coincide with a point   on the tesselated
    % outer surface layer grid, the sensor location will be quantized to the
    % closest outer surface tesselation point. Dipole generator(s) are assumed to be
    % interior  to the innermost "core" layer
    %
    % INPUTS (Required)
    %       R_eeg  : EEG sensor locations on scalp (meters)         (Meeg x 3)
    %                  (insert dummy value if mode = 2)
    %       R_meg  : MEG sensor locations on scalp (meters)         (Mmeg x 3)
    %                  (insert dummy value if mode = 1)
    %       O_meg  : MEG sensor orientations (unit-vector)          (Mmeg x 3)
    %                  (insert dummy value if mode = 1)
    %    vertices  : Tessalated Surface layers from INNERMOST TO OUTERMOST.
    %                Packed as Cell Array where index 1 is innermost surface  (1 x NL)-Cell Array
    %                Individual entries are position row vectors in PCS
    %       faces  : Tesselated Surface Node connections for each triangle    (1 x NL)-Cell Array
    %                Packed as Cell Array where index 1 is innermost surface
    %                Individual entries are indexes to "vertices"
    %       sigma  : conductivity from INNERMOST to OUTERMOST        (NL x 1)
    %       mode   : 1,  compute EEG only
    %                2,  compute MEG only
    %                3,  compute both EEG and MEG
    %   basis_opt  : 0, constant basis
    %                1, linear
    %   test_opt   : 0, collocation,
    %                1, Galerkin
    %     ISA      : 0,  Inhibit Isolated Skull Approach
    %              : 1,  Enable Isolated Skull Approach
    %
    %
    %    fn_eeg    : EEG transfer matrix filename                   Char String
    %                   - EEG Transfer matrices are computed and then saved
    %                     in this "*.mat" file
    %
    %    fn_meg    : MEG transfer matrix filename                   Char String
    %                   - MEG Transfer matrices are computed and then saved
    %                     in this "*.mat" file
    %
    %         WHERE: M = # of sensors; P = # of dipoles; NL = # of sphere layers
    %               Nn = # of Surface Tesselation Nodes; Nt = # of Surf Tess Triangles
    %
    % INPUTS (Optional):
    %    NVertMax : Maximum Number of vertices for Surface Tesselation  (for a single layer)
    %                   - Surface is re-tesselated if # points exceeds max
    %                   - Default Maximum Number of Triangle Faces is 2292    (1 x NL)
    %
    %    Verbose : Toggles Verbose mode on/off (on = 1);
    %
    % OUTPUTS:
    %    nfv    : a surface tessellation structure with fields .vertices and
    %             .faces of all envelopes used in BEM computation, once
    %             downsampling, separation and embedding check have been
    %             completed.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %<autobegin> ---------------------- 14-Jun-2004 17:09:48 -----------------------
    % --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
    %
    % CATEGORY: Forward Modeling
    %
    % Alphabetical list of external functions (non-Matlab):
    %   publictoolbox\othertools\inpolyhd.m
    %   toolbox\bst_message_window.m
    %   toolbox\crossprod.m
    %   toolbox\dotprod.m
    %   toolbox\get_user_directory.m
    %   toolbox\hrow_linear.m
    %   toolbox\makeuswait.m
    %   toolbox\ray_intersect.m
    %   toolbox\rownorm.m
    %   toolbox\solid_angle2.m
    %   toolbox\tessellation_outwards.m
    %   toolbox\trans_matrix.m
    %   toolbox\view_surface.m
    %
    % Subfunctions in this file, in order of occurrence in file:
    %   [H,A,G] = bem_linear(mode,choice,cdvdiff,cdvsum,geometry,...
    %   [H,A] = bem_constant(mode,choice,cdvdiff,cdvsum,r1,r2, ...
    %   I = nit_gq(F,N,Idim,r1,r2,r3,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12)
    %   nfv = make_surfaces_coincide(fv,Verbose)
    %   fv = make_surfaces_embedded(fv,tol,Verbose)
    %
    % At Check-in: $Author: Mosher $  $Revision: 23 $  $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:48 -----------------------

    % Older comments about outputs from original bem_xfer version (J. Ermer)

    %(Optional: Complete set of output parameters are stored in
    %                     user specified "*.mat" files )
    %    Te_ISA    : EEG BEM transfer matrix (Isolated Skull Approach)
    %    Te        : EEG BEM transfer matrix
    %    Tm_ISA    : MEG BEM transfer matrix (Isolated Skull Approach)
    %    Tm        : MEG BEM transfer matrix
    %
    % ----------------------------- Script History ---------------------------------
    % External Functions and Files:
    %    bem1.m: Computes BEM transfer matrices (James Chang 1995)
    %    rownorm.m; dotprod.m: USC/LANL MEG/EEG Toolbox
    %    kernel.m: Computes Forward Gain Matrices using BEM transfer matrices (James Chang 1995)
    %
    %     - John J. Ermer 6/3/99
    %        -  2/21/00: Modified to compute EEG/MEG transfer matrices only (J. Ermer)
    %        -  2/25/00: Integrate James Chang Routine directly into code (J. Ermer)
    %        -  2/29/00: Added capability for independent # of vertices for each surface layer
    %        -  3/01/00: Modified Inputs to be passed in Cell Array Format (J. Ermer)
    %        -  3/16/00: Added minimum memory requirement model (J. Ermer)
    %        -  4/23/00: Added logic to allow for interpolation across triangle surface
    %                    for Linear BEM (J. Ermer)
    %
    % Date of Modifiation: March, 2002
    % March 25, 2002 : SB: added control on Verbose mode
    % SB     May-2004  Extensive upgrades, checks
    % SB  28-May-2004  Added ray_tracing for coincident vertices, embedding, upgraded
    %                  detection of memory requirements to account for
    %                  constant/linear cases, altered reducepatch calls.
    % JCM 28-May-2004  Editing, cleaning comments, bug chasing.
    % ----------------------------- Script History ---------------------------------

    %
    %%% THIS PART DECLARES PROGRAM CONSTANTS AND PARAMETERS %%%
    %

    if Verbose
        bst_message_window(sprintf('Computing BEM Transfer Matrices. . .'))
    end

    Meeg = size(R_eeg,1);   % Number of EEG Sensors
    Mmeg = size(R_meg,1);   % Number of MEG Sensors
    NL = size(vertices,2);  % Number of Surface Layers

    % MaxMemoryConstraint
    MaxMemoryConstraint = 500; %MB

    % Now depending on the basis, test how many vertices per envelope we need
    memflag = 0;
    switch basis_opt
    case 0 % Constant basis
        MemLoad = (NL*(2*NVertMax-4))^2*8/1e6; % in MB / Number of Faces is the parameter
        if MemLoad > MaxMemoryConstraint | ~exist('NVertMax','var')
            NVertMax = round(.5*( (1/NL)*sqrt(1e6*MaxMemoryConstraint/8) + 4));
            memflag = 1;
        end

    case 1 % Linear
        MemLoad = (NL*NVertMax)^2*8/1e6; % in MB
        if MemLoad > MaxMemoryConstraint | ~exist('NVertMax','var')
            NVertMax = round((1/NL)*sqrt(1e6*MaxMemoryConstraint/8));
            memflag = 1;
        end
    end

    if memflag
        bst_message_window('wrap',...
            {...
                'Warning: possible memory limitation detected',...
                sprintf('Reducing number of vertices to approx. %d',NvertMax)...
            })
    end

    %
    %%% THIS PART CHECKS INPUT PARAMETERS FOR THEIR DIMENSION AND VALIDITY %%%
    %
    %
    if (size(R_eeg,1)>=1)&(size(R_eeg,2)~=3)
        error('EEG sensor location(s) must have three columns!!!')
    end
    %
    if (size(R_meg,1)>=1)&(size(R_meg,2)~=3)
        error('MEG sensor location(s) must have three columns!!!')
        if max(abs(rownorm(O_meg)-ones(Mmeg,1))>1.0e-3)
            error('MEG Sensor Orientations are not Unity!!!')
        end
        if size(O_meg,1)~=size(R_meg,1)
            error('Number of MEG Sensors must equal number of MEG Sensor Orientations!!!')
        end
    end
    %
    if (Mmeg>=1)&(size(O_meg,2)~=3)
        error('MEG sensor orientation(s) must have three columns!!!')
    end
    %
    %%%% THIS PART CHECKS THE DIMESION OF PARAMETERS ASSOCIATED WITH EACH SURFACE LAYER %%%%
    %
    for k=1:NL
        if size(vertices{k},2)~=3
            vertices{k} = vertices{k}';
        end

        if size(faces{k},2)~=3
            faces{k} = faces{k}';
        end

        if size(faces)~=size(vertices)
            error('Variables "vertices" and "faces" must have same number of layers (NL)!!!')
        end

        Nv(k) = size(vertices{k},1);
        % Number of Surface Node (tesselation) Points (for each surface)

        Nf(k) = size(faces{k},1);
        % Number of Surface Triangles (for each surface

    end
    %
    %%% THIS PART DOWNSAMPLES THE NUMBER OF TESSELATION NODES IF THE NUMBER OF POINTS EXCEEDS A %%%%%%%%%
    %%% USER SPECIFIED -OR- DEFAULT MAXIMUM %%%%%%
    %
    % if ~exist('NVertMax','var') | isempty(NVertMax),
    %     NVertMax = NVertMax_default;
    % end
    %
    for k=1:NL
        fv(k).faces = double(faces{k});
        fv(k).vertices = double(vertices{k});

        if Nf(k) > NVertMax
            if Verbose & k==1
                bst_message_window(sprintf('Reducing all %d tessellated layers to %d vertices. . .',NL,NVertMax));
            end
            % Reduce tesselation to required maximum number of vertices
            % (call to reducepatch/qslim necessitates specification of the number of faces though)
            [fv(k)] = reducepatch(fv(k),2*NVertMax-4);
            %[fv(k)] = qslim(fv(k),'-t',2*NVertMax-4);

            % Check for triangle orientation (need to point outwards)
            [fv(k), status] = tessellation_outwards(fv(k));

            if ~status & Verbose
                bst_message_window(sprintf('Surface #%d had wrong triangle orientation; changed all.',k))
            end

            Nf(k) = size(fv(k).faces,1);                  % Reduced number of triangular faces
            Nv(k) = size(fv(k).vertices,1);                   % Reduced number of vertices
            if Verbose & k==NL
                bst_message_window(sprintf('Reducing all %d tessellated layers to %d vertices -> DONE',NL,NVertMax));
            end
        end
    end

    if Verbose
        % Now realign surface nodes (SB Feb 18, 2003)
        bst_message_window('wrap',...
            'Aligning surface nodes for coincident tessellated envelopes. . .');
    end

    nfv = make_surfaces_coincide(fv,Verbose);
    %nfv = make_collinear_surface2(fv,0,Verbose);

    tol = .003; %m % Minimum distance between surfaces
    if Verbose
        % Now realign surface nodes (SB Feb 18, 2003)
        bst_message_window('wrap',...
            sprintf('Checking for minimum separation of %3.1f mm between envelopes. . .',1000*tol));
    end

    nfv = make_surfaces_embedded(nfv,tol,Verbose);

    if Verbose
        % Now realign surface nodes (SB Feb 18, 2003)
        bst_message_window('wrap',...
            '-> DONE');
    end

    % Final adjustment
    nfv = make_surfaces_coincide(nfv,Verbose);
    if nargout > 0
        varargout{1} = nfv;
    end

    if Verbose

        bst_message_window('wrap',...
            'Aligning surface nodes for collinear tessellated envelopes -> DONE');

        % Visualization of surfaces + grid points
        for k = 1:length(fv)
            [hf,hs(k),hl] = view_surface('Original Surfaces',fv(k).faces,fv(k).vertices);
            [nhf,nhs(k),nhl] = view_surface('Realigned Surfaces2',nfv(k).faces,nfv(k).vertices);
        end
        rotate3d on

        set(hs(1),'FaceAlpha',.8,'edgealpha',1,'edgecolor','r','facecolor','r')
        set(hs(2),'FaceAlpha',.2,'edgealpha',1,'edgecolor','g','facecolor','g')
        set(hs(3),'FaceAlpha',.1,'edgealpha',1,'edgecolor','b','facecolor','b')
        set(nhs(1),'FaceAlpha',.8,'edgealpha',1,'edgecolor','r','facecolor','r')
        set(nhs(2),'FaceAlpha',.2,'edgealpha',1,'edgecolor','g','facecolor','g')
        set(nhs(3),'FaceAlpha',.1,'edgealpha',1,'edgecolor','b','facecolor','b')

    end

    for k = 1:length(nfv)
        faces{k} = nfv(k).faces;
        vertices{k} = nfv(k).vertices;
    end
    clear nfv fv

    %
    %%%% THIS PART DETERMINES ESTIMATED MEMORY REQUIREMENTS FOR COMPUTING THE TRANSFER MATRICES %%%%%%%
    %
    if basis_opt==0  % Constant Case
        Nftot = sum(Nf);
        mem_min = (3*Nftot*Nftot ....        % H,L,U Matrices
            + 3*Mmeg*Nftot ....         % A Matrix
            + 11*Nftot*3 ....           % Misc Workspace Matrices
            + 12*Nftot*1)*8*1.0e-6;     % Misc Workspace Matrices
    elseif basis_opt==1  % Linear Case
        Nvtot = sum(Nv);
        mem_min = (3*Nvtot*Nvtot ....        % H,L,U Matrices
            + 3*Mmeg*Nvtot ....         % A Matrix
            + 11*Nvtot*3 ....           % Misc Workspace Matrices
            + 12*Nvtot*1)*8*1.0e-6;     % Misc Workspace Matrices
    end
    if Verbose
        bst_message_window('wrap',{...
                sprintf('Updated estimated virtual memory load (including misc. worskspace matrices) : %.3f MB.',mem_min),...
            })
    end

    %
    %%%% This part formats parameters for pre-existing routines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% NOTE: GEOMETRY AND NODES REFORMATTED TO OUTERMOST->INNERMOST FORMAT %%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    geometry = zeros(sum(Nf),3);   % Triangle Interconnection Matrix for each layer
    % specifying nodes for each triangle
    nodes = zeros(sum(Nv),4);      % Nodes for each layer; Col 4 = layer # (Outermost=1!!!)
    %
    try
        for k=1:NL  % Reverse order from outermost to innermost !!!
            geometry((sum(Nf(NL-k+2:NL))+1):(sum(Nf(NL-k+1:NL))),1:3) = faces{NL-k+1} + repmat(sum(Nv(NL-k+2:NL)),Nf(NL-k+1),3);
            nodes((sum(Nv(NL-k+2:NL))+1):(sum(Nv(NL-k+1:NL))),1:4) = [vertices{NL-k+1} k*ones(Nv(NL-k+1),1)];
        end
    catch
        errordlg('Surfaces have probably too few vertices. Try to increase the number of vertices per tessellated envelope.', mfilename)
        return
    end

    %
    r1 = nodes(geometry(:,1),1:3);   % All Layers Vertex 1
    r2 = nodes(geometry(:,2),1:3);   % All Layers Vertex 2
    r3 = nodes(geometry(:,3),1:3);   % All Layers Vertex 3
    ctrd = (r1+r2+r3)/3;         % Centroid of each triangle (all layers)
    %
    % compute the normal vectors
    tri_xp=(cross((r2-r1)',(r3-r1)'))';  % Triangle cross product
    area2=rownorm(tri_xp);            % twice area of triangle
    N= tri_xp./ [area2,area2,area2];  % normalization
    %
    %%%% FOR THE COLLOCATION CASE ONLY, THIS PART FINDS THE TRIANGLE ON THE OUTERMOST LAYER %%%%%%%%%%%
    %%%% WHOSE CENTROID IS CLOSEST TO EACH EEG SENSOR. THE SENSOR LOCATION IS THEN QUANTIZED %%%%%%%%%%
    %%%% TO THE LOCATION OF THE CLOSEST TRIANGLE CENTROID %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    P = zeros(Meeg,1);
    P_sens = zeros(Meeg,3);
    P_wts = [];
    if Meeg >= 1
        for m=1:Meeg
            [temp indx] = min(rownorm(repmat(R_eeg(m,:),Nf(NL),1)- ctrd(1:Nf(NL),:)));
            P(m) = indx;                        % Index of Tesselation Triangle in "geometry"
            P_sens(m,:) = ctrd(indx,:);
            R_eeg_dist(m,:) = temp;
            R_eeg_quant(m,:) = ctrd(indx,:);
        end
    else
        P=[];
    end
    %
    %%%% FOR THE LINEAR CASE ONLY, THE THREE VERTICES ASSOCIATED WITH THE TRIANGLE WHOS CENTROID %%%%%%
    %%%% (COMPUTED ABOVE)IS CLOSEST TO THE SPECIFIED POSITION ARE FOUND. DURING KERNEL OPERATIONS, %%%%
    %%%% THE FINAL SOLUTION IS INTERPOLATED BETWEEN THE 3 VERTICES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    if basis_opt==1  % Linear Case
        P = zeros(3*Meeg,1);       % Node index list
        P_sens = zeros(3*Meeg,3);    % Node Position List
        P_wts = zeros(Meeg,3);     % Node Weight List (for sensor position lying interior to triangle)
        if Meeg >= 1
            %
            for m=1:Meeg  % Loop thourgh sensors
                %
                %%%% This part finds closest node to sensor position and those triangles containing node %%%%
                %
                [temp indx] = min(rownorm(repmat(R_eeg(m,:),Nv(NL),1)- nodes(1:Nv(NL),1:3)));
                %
                Pt =[];
                for  i=1:3,                            % find triangles adjacent to node p
                    t =find(geometry(:,i)==indx);
                    Pt= [Pt; t];  % Append to adjacent triangle list
                end
                %
                %%%% This part finds the projection of sensor point onto nearest triangle surface %%%%
                for i=1:size(Pt,1)  % check each triangle to find the one which contains projected sensor
                    Atemp = [r2(Pt(i,1),:)-r1(Pt(i,1),:); r3(Pt(i,1),:)-r1(Pt(i,1),:)]'; % triangle basis
                    btemp = R_eeg(m,:)' - r1(Pt(i,1),:)';    % Sensor point
                    xtemp = Atemp\btemp;   % Solve for linear combination
                    pnt = r1(Pt(i,1),:) + (Atemp*xtemp)';   % Projection of sensor point onto triangle surface
                    %%%% This part computes the weights for each vertices; If point outside triangle %%%%%%%%%%%%%%%%%%
                    %%%% nearest vertice node is given full weight %%%%
                    %
                    P_wts(m,3) = N(Pt(i,1),:)*crossprod(r2(Pt(i,1),:)-r1(Pt(i,1),:),pnt-r1(Pt(i,1),:))'/area2(Pt(i,1),:); %s
                    P_wts(m,2) = N(Pt(i,1),:)*crossprod(pnt-r1(Pt(i,1),:),r3(Pt(i,1),:)-r1(Pt(i,1),:))'/area2(Pt(i,1),:); %t
                    P_wts(m,1) = 1 - P_wts(m,2) - P_wts(m,3);
                    %
                    P(3*m-2:3*m,1) = geometry(Pt(i),1:3)';  % Node List (Triangle Vertices)
                    P_sens(3*m-2:3*m,:) = nodes(geometry(Pt(i),1:3)',1:3);
                    %
                    if (min(P_wts(m,:)')>=0)&(max(P_wts(m,:)')<=1) % Stop if pnt falls interior to triangle
                        break
                    end
                    %
                end  % Triangle List Loop
                %
            end  % Sensor Loop
        else
            P=[];
        end
        %
    end
    %
    %%%% IF NO MEG SENSORS HAVE BEEN SPECIFIED, THIS PART GENERATES A DUMMY SENSOR TO AVOID ERROR %%%%
    %%%% STATEMENTS IN EXISTING CODE (TEMPORARY MEASURE)
    %
    if Mmeg >= 1
        S = R_meg;
    else
        S = [];
    end
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% THIS PART COMPUTES THE BEM TRANSFER MATRICES BASED ON JAMES CHANG'S ORIGINAL ROUTINE %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %

    %
    t0 = clock;
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    %%%% THIS PART FORMATS CONDUCTIVITY PARAMETERS %%%%
    %
    len = length(sigma);
    cdvity = [sigma,0];
    cdv =zeros(len,2);

    for i=len:-1:1,
        cdv(len+1-i,:)=[cdvity(i),cdvity(i+1)];
    end

    ss=cdv(:,1)+cdv(:,2);   % conductivity sum no_shell x 1
    d=cdv(:,1)-cdv(:,2);   % conductivity difference no_shell x 1
    %
    %%%% THIS PART COMPUTES THE NORMAL VECTORS FOR EACH TRIANGLE %%%%
    %
    % change directions to be outward
    dp=(dot(ctrd',N'))';
    dex=find(dp < 0);
    N(dex,:)= -N(dex,:);
    temp=r2(dex,:);      %  r1,r2,r3 is CCW if looked outside the head
    r2(dex,:)=r3(dex,:);
    r3(dex,:)=temp;
    %
    temp=geometry(dex,2);      %  CCW
    geometry(dex,2)=geometry(dex,3);
    geometry(dex,3)=temp;
    no_elt = size(geometry,1);
    no_node = size(nodes,1);
    %
    clear temp dp dex
    %
    if basis_opt == 0,
        whichsurf = nodes(geometry(:,1),4);
    elseif basis_opt == 1 | basis_opt==2,
        whichsurf = nodes(:,4);
    end
    %
    cdvsum=ss(whichsurf);
    cdvdiff=d(whichsurf);
    no_surf =max(whichsurf);
    DIM =zeros(no_surf,1);
    for i=1: no_surf
        DIM(i) = length(find( whichsurf ==i));
    end
    %
    %%%% THIS PART COMPUTES THE H AND A MATRICES %%%%
    %
    H = []; A = [];
    if mode ==1, S=[]; end
    if basis_opt == 0         % CONSTANT COLLOCATION CASE
        if Verbose
            bst_message_window('-> BEM with constant basis') % H : no of triangle x no of triangle
        end

        eegsensor = ctrd(P,1:3);     % EEG SENSOR LOCATION SET TO TRIANGLE CENTROID
        [H,A] = bem_constant(mode,test_opt,cdvdiff,cdvsum,r1,r2,r3,ctrd,area2/2,N,S);
        %  [H,A] = bem_constant2(mode,test_opt,cdvdiff,cdvsum,r1,r2,r3,ctrd,area2/2,N,S);
        I = speye(no_elt);

    elseif basis_opt == 1    % LINEAR COLLOCATION CASE
        if Verbose
            bst_message_window('-> BEM with linear basis')
        end
        % H: no of nodes x no of nodes
        eegsensor = nodes(P,1:3);   % EEG SENSOR LOCATION SET TO CLOSEST NODE
        [H,A,I] = bem_linear(mode,test_opt,cdvdiff,cdvsum,geometry,...
            nodes,r1,r2,r3,N,area2,S);
    end
    %
    %%%% This part does other functions
    %
    PP =  sparse(length(P),sum(DIM)); % a sparse zero matrix
    for i=1:length(P);
        PP(i,P(i))=1;
    end
    %
    if mode ==1, S=[]; end;

    if Verbose
        bst_message_window('Completing LU decomposition. . .'), end

    [Te_ISA,Te,Tm_ISA,Tm]= trans_matrix(mode,test_opt,ISA,DIM,cdv(no_surf,:),I,H,PP,A,size(S,1));

    if Verbose
        bst_message_window('Completing LU decomposition -> DONE'), end

    %
    Te_ISA = full(Te_ISA);
    Te = full(Te);
    Tm_ISA = full(Tm_ISA);
    Tm = full(Tm);
    %
    basisname = ['constant BEM';...
            'linear   BEM'];
    testname = ['collocation  ';...
            'galerkin     '];
    basis = basisname(basis_opt+1,:);
    test = testname(test_opt+1,:);
    %
    telap_xfer = etime(clock,t0);
    if Verbose
        bst_message_window(sprintf('Computed in %3.1f sec.',telap_xfer))
    end

    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% THIS PART STORES FINAL TRANSFER MATRIX RESULTS AND PARAMETERS IN A ".MAT" FILE %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    eegsensor = [P(:),P_sens];
    if (mode==1)|(mode==3)
        modality = 1;
        try
            eval(['save ', fn_eeg,[' basis_opt test_opt basis test geometry nodes cdv' ...
                        ' Te Te_ISA eegsensor modality R_eeg P_sens P_wts']]);
        catch
            Users = get_user_directory;
            cd(Users.STUDIES)
            eval(['save ', fn_eeg,[' basis_opt test_opt basis test geometry nodes cdv' ...
                        ' Te Te_ISA eegsensor modality R_eeg P_sens P_wts']]);
        end

        if Verbose
            bst_message_window(sprintf('EEG Transfer Matrix data stored in file %s.',fn_eeg))
        end

        %
    end

    if  (mode==2)|(mode==3)
        modality = 2;
        try
            eval(['save ', fn_meg,[' basis_opt test_opt basis test geometry nodes cdv' ...
                        ' Tm Tm_ISA S modality R_meg O_meg']])
        catch
            Users = get_user_directory;
            cd(Users.STUDIES)
            eval(['save ', fn_meg,[' basis_opt test_opt basis test geometry nodes cdv' ...
                        ' Tm Tm_ISA S modality R_meg O_meg']])
        end

        if Verbose
            bst_message_window(sprintf('MEG Transfer Matrix data stored in file %s.',fn_meg))
        end

        %
    end
    %
    if Verbose
        bst_message_window(sprintf('Computing BEM Transfer Matrices -> DONE'))
    end

    % ----------------------------------------------------------------------------------------------------------------

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    %
    % SUBFUNCTIONS

    % ----------------------------------------------------------------------------------------------------------------

    function [H,A,G] = bem_linear(mode,choice,cdvdiff,cdvsum,geometry,...
        nodes,r1,r2,r3,N,area2,S)
    %bem_linear (overwrite succinct one line summary here)
    % function [H,A,G] = bem_linear(mode,choice,cdvdiff,cdvsum,geometry,...
    %                             nodes,r1,r2,r3,N,area2,S)
    %
    % r1(i,:), r2(i,:), r3(i,:) form three vertices of triangle i and ther are
    % CCW as seen from outside
    % References:
    % (1) de Munck,IEEE Trans. BME, pp 986-990,1992
    % (2) Schlitt,et al, IEEE Trans. BME, pp 52-58,1995
    %
    % %%% (03/12/00)-Added Case where MEG is not computed to prevent
    %                Matlab error from being flagged (John Ermer)
    %

    %<autobegin> -------- 20-Nov-2002 14:04:25 ------------------------------
    % ---- Automatically Generated Comments Block using auto_comments -----------
    %
    % Alphabetical list of external functions (non-Matlab):
    %   toolbox\crossprod.m
    %   toolbox\dotprod.m
    %   toolbox\hrow_linear.m
    %   toolbox\nit_gq.m
    %   toolbox\rownorm.m
    %<autoend> ---------- 20-Nov-2002 14:04:25 ------------------------------

    no_node = size(nodes,1);
    no_elt=size(geometry,1);
    no_eltx2 =2*no_elt;

    PP =  sparse(no_node,3*no_elt);
    for i=1:no_elt,
        a=  1/area2(i);
        PP(geometry(i,1),i)=a;
        PP(geometry(i,2),i+no_elt)=a;
        PP(geometry(i,3),i+no_eltx2)=a;
    end
    % work on 'geometry' matrix.
    % three vertices of the ith row contribute to the ith triangle i.

    r21 = r2-r1;
    r32 = r3-r2;
    r13 = r1-r3;
    yset = [r21,r32,r13];
    yset_mag = [rownorm(r21), rownorm(r32),rownorm(r13)];

    dx = [ 2 3 1 2];
    dx1= [1 2 3 1];

    H = zeros(no_node);
    %disp('Computing H .... ')

    if choice == 0, % collocation
        G = speye(no_node);
        %  t=  clock;
        for  i=1:no_node
            v =1:no_elt;
            x =find(geometry(:,1)==i | geometry(:,2)==i | geometry(:,3)==i);
            % excluding all triangles having node i as a vertex as the integral
            % across these trianglesa are either zeros or singular
            v(x)=[];
            H(i,:) = hrow_linear(nodes(i,1:3),PP,v,yset,yset_mag,geometry,nodes,cdvdiff,N,area2)/(2*pi*cdvsum(i));

        end

        for i=1:no_node,    % diagonal terms
            H(i,i)= 1-sum(H(i,:));
        end

    else
        G = sparse(no_node,no_node);
        rnxrn1 = [crossprod(r2,r3),crossprod(r3,r1),crossprod(r1,r2)];
        DET = dotprod(r1,rnxrn1(:,1:3));
        for p=1:no_node,
            %    tt=  clock;
            %sprintf('node# = %5.0f',p)
            %       profile on -detail builtin
            ptri =[];
            for  i=1:3,   % find triangles adjacent to node p
                t =find(geometry(:,i)==p);
                ptri= [ptri; t,i*ones(size(t))];
            end

            for i=1:size(ptri,1)  % compute the Gram matrix
                t= ptri(i,1);
                G(p,geometry(t,:)) = G(p,geometry(t,:))+ nit_gq('kl_gram',4,3,r1(t,:),r2(t,:),r3(t,:),rnxrn1(t,:),ptri(i,2))/DET(t)^2;
            end

            for i=1:size(ptri,1),
                t = ptri(i,1);
                v = 1:no_elt;
                v(t)=[];
                H(p,:) = H(p,:)+nit_gq('kl_galerkin',3,no_node,r1(t,:),r2(t,:),...
                    r3(t,:),PP,ptri(i,2),rnxrn1(t,:),v,yset,yset_mag,...
                    geometry,nodes,cdvdiff,N,area2)/DET(t);
            end

            H(p,:) = H(p,:)/(2*pi*cdvsum(p));    % cdvdiff is treated in hrow_linea
            %profile report
        end
    end    % end if

    %
    % ----    compute the A matrix for MEG exactly   -----------
    % ----    see Ferguson's paper for details ------------------
    %
    if mode==2 | mode ==3,
        %disp('Computing the matrix A for MEG (omitting u0/4pi) .... ')
        % t=clock;
        if length(S)==3,
            no_megsensor=1;     % only one meg sensor
            S=S(:)';            % into a row vector
        else  no_megsensor=size(S,1);
        end

        A=zeros(3*no_megsensor,no_node);
        cols = [0,no_elt,no_elt*2,no_elt*3];
        for i=1:no_megsensor,
            y1 = [r1(:,1)-S(i,1),r1(:,2)-S(i,2),r1(:,3)-S(i,3)];
            y2 = [r2(:,1)-S(i,1),r2(:,2)-S(i,2),r2(:,3)-S(i,3)];
            y3 = [r3(:,1)-S(i,1),r3(:,2)-S(i,2),r3(:,3)-S(i,3)];
            mag1=rownorm(y1);
            mag2=rownorm(y2);
            mag3=rownorm(y3);
            sangle = 2*atan2(dotprod(y1,crossprod(y2,y3)),mag1.*mag2.*mag3+...
                mag1.*dotprod(y2,y3) + mag2.*dotprod(y1,y3) +...
                mag3.*dotprod(y1,y2) );

            c = dotprod(y1,N);
            c = N.*[c,c,c];
            a = [S(i,1)+c(:,1),S(i,2)+c(:,2),S(i,3)+c(:,3)];
            c1 = r1-a;
            c2 = r2-a;
            c3 = r3-a;

            tmp = zeros(no_elt,3);
            tmpsum =zeros(no_elt,1);

            for p=1:3,
                i1 =dx1(p+1);
                i2 = p;
                yp1 = eval(['y',num2str(i1)]);
                yp = eval(['y',num2str(i2)]);
                cp1 = eval(['c',num2str(i1)]);
                cp = eval(['c',num2str(i2)]);

                yp1_yp = yset(:,3*p-2:3*p);
                mag_yp1_yp = yset_mag(:,p);

                w = (mag_yp1_yp.*rownorm(yp1) + dotprod(yp1,yp1_yp))./...
                    (mag_yp1_yp.*rownorm(yp) + dotprod(yp,yp1_yp));
                gammap =log(w)./mag_yp1_yp;
                tmp =tmp+crossprod(cp,cp1).*[gammap,gammap,gammap];
            end
            % twice area is absorbed into PP
            tmpsum =( dotprod(N,tmp)-dotprod(N,c).*sangle);
            rows = 3*i-2:3*i;
            for n=1:3,
                col= geometry(:,n);
                yj_yk = -yset(:,3*dx(n)-2:3*dx(n));

                a = cdvdiff(col).*tmpsum;
                a = (yj_yk .*[a,a,a]);

                A(rows,:)= A(rows,:)+(PP(:,cols(n)+1:cols(n+1))*a)';
            end

        end
    else
        A = [];
    end

    % ----------------------------------------------------------------------------------------------------------------

    function [H,A] = bem_constant(mode,choice,cdvdiff,cdvsum,r1,r2, ...
        r3,ctrd,area,N,S)
    %bem_constant Compute the geometry matrix H for EEG & MEG by centroid approxiamtion
    % function [H,A] = bem_constant(mode,choice,cdvdiff,cdvsum,r1,r2, ...
    %                  r3,ctrd,area,N,S)
    % - compute the geometry matrix H for EEG & MEG by centroid approximation of numerical integration
    % - compute the matrix A analytic
    % See:
    % Hamalainen, Sarvas, IEEE Trans. BME, pp.165-171, 1989
    % de Munck,IEEE Trans.BME,pp986-990,1992 (for integration of the MEG kernel)
    %
    % (03/12/00)-Added Case where MEG is not computed to prevent
    %                Matlab error from being flagged (John Ermer)

    %<autobegin> -------- 20-Nov-2002 14:04:16 ------------------------------
    % ---- Automatically Generated Comments Block using auto_comments -----------
    %
    % Alphabetical list of external functions (non-Matlab):
    %   toolbox\dotprod.m
    %   toolbox\nit_gq.m
    %   toolbox\rownorm.m
    %   toolbox\solid_angle2.m
    %<autoend> ---------- 20-Nov-2002 14:04:16 ------------------------------

    %TH = max(rownorm(r1)); % threshold used in Galerkin method, see below
    %tic
    no_elt = size(r1,1);

    %%disp(' Computing  H ......')
    H=zeros(no_elt);

    for i=1:no_elt
        %sprintf('Node #%4f...',i)
        x=1:no_elt;
        x(i)=[];     % excluding the diagonal terms
        sangle = zeros(no_elt-1,1);
        if choice ==0,    % ctrd approxmation
            sangle = solid_angle2(ctrd(i,:),r1(x,:),r2(x,:),r3(x,:));

        else  % Galerkin
            sangle =nit_gq('solid_angle2',3,1,r1(i,:),r2(i,:),r3(i,:),r1(x,:),r2(x,:),r3(x,:))/area(i);

        end      % end if
        H(i,x) = ( cdvdiff(x).* sangle )'/( 2*pi*cdvsum(i)) ;
    end;
    %telap_h1 = toc
    %
    %  compute the A matrix for MEG exactly
    %  see de Munck's paper for details
    %
    if mode==2 | mode ==3 , % MEG and Fusion case
        %disp('Computing the matrix A for MEG(omitting u0/4pi) .... ')
        % t = clock;
        if length(S)==3,
            no_megsensor=1;     % only one meg sensor
            S=S(:)';            % into a row vector
        else  no_megsensor=size(S,1);
        end

        A=zeros(3*no_megsensor,no_elt);
        r1_r3= r1-r3;
        r3_r2= r3-r2;
        r2_r1= r2-r1;
        mag_r2_r1=rownorm(r2_r1);
        mag_r3_r2=rownorm(r3_r2);
        mag_r1_r3=rownorm(r1_r3);

        for i=1:no_megsensor,
            S_r1=[S(i,1)-r1(:,1),S(i,2)-r1(:,2),S(i,3)-r1(:,3)];
            S_r2=[S(i,1)-r2(:,1),S(i,2)-r2(:,2),S(i,3)-r2(:,3)];
            S_r3=[S(i,1)-r3(:,1),S(i,2)-r3(:,2),S(i,3)-r3(:,3)];
            mag_S_r1=rownorm(S_r1);
            mag_S_r2=rownorm(S_r2);
            mag_S_r3=rownorm(S_r3);

            gamma1 = -log(( mag_r2_r1.*mag_S_r2 - dotprod(r2_r1,S_r2)) ...
                ./( mag_r2_r1.*mag_S_r1-dotprod(r2_r1,S_r1) )) ./mag_r2_r1;

            gamma2 = -log(( mag_r3_r2.*mag_S_r3 - dotprod(r3_r2,S_r3)) ...
                ./( mag_r3_r2.*mag_S_r2-dotprod(r3_r2,S_r2) )) ./mag_r3_r2;

            gamma3 = -log(( mag_r1_r3.*mag_S_r1 - dotprod(r1_r3,S_r1)) ...
                ./( mag_r1_r3.*mag_S_r3-dotprod(r1_r3,S_r3) )) ./mag_r1_r3;

            % omit u0/4pi
            A(3*i-2:3*i,:)= (([gamma1,gamma1,gamma1].*r2_r1+...
                [gamma2,gamma2,gamma2].*r3_r2+...
                [gamma3,gamma3,gamma3].*r1_r3) .*[cdvdiff,cdvdiff,cdvdiff])';
        end
    else % EEG case
        A = [];
    end
    %keyboard

    % ----------------------------------------------------------------------------------------------------------------

    function I = nit_gq(F,N,Idim,r1,r2,r3,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12)
    %NIT_GQ (F,N,Idim,r1,r2,r3,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12)
    % function I = nit_gq(F,N,Idim,r1,r2,r3,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12)
    % This function integrates a scalar/vector(dimension indicated by Idim)
    % function over triangles with vertices r1,r2,and r3, using
    % the Gauss Quadrature,i.e.,
    %    /
    %    | F(r,P1,P2,P3,....)da    ,
    %    /
    % where  r=(x,y,z) is a point on the triangle, and P1,P2,.... are parameters
    % defined in F.
    %
    % Input:
    % F     :  A string containing the function to be integrated. It should take r,P1,P2,P3,...
    %          as input, where r is no_triangles by 3,and produces no_triangles by
    %          Idim funtion valus.
    % N     :  use N point Gauss Quadrature          scalar
    % r1, r2,r3 :  specify the vertices of triangles over which F is integrated;
    %              each is no_triangles x 3
    % P1,P2..  : parameters defined in F.
    %
    % Output:
    % I : integral                            no_triangles(=size(r1,1)) x Idim
    %

    %<autobegin> -------- 20-Nov-2002 14:07:22 ------------------------------
    % ---- Automatically Generated Comments Block using auto_comments -----------
    %
    % Alphabetical list of external functions (non-Matlab):
    %   toolbox\crossprod.m
    %   toolbox\rownorm.m
    %
    %<autoend> ---------- 20-Nov-2002 14:07:22 ------------------------------

    % Author: CCH
    %

    if N < 2 | N > 10,
        error(' The GQ table of your choice is not built. Sorry');
    end

    %-----------------------------------------------------------
    % W: weight table, X : abscissa table
    X2 = 0.5773503;
    W2 = 1.0;

    X3 = [0.0,       0.7745967];
    W3 = [0.4444445, 0.5555556];

    X4 = [0.3399810, 0.8611363];
    W4 = [0.6521452, 0.3478548];

    X5 = [0.0,       0.5384693, 0.9061798];
    W5 = [0.2844445, 0.4786287, 0.2369269];

    X6 = [0.2386192, 0.6612094, 0.9324695];
    W6 = [0.4679139, 0.3607616, 0.1713245];

    X7 = [0.0,       0.4058452, 0.7415312, 0.9491079];
    W7 = [0.2089796, 0.3818301, 0.2797054, 0.1294850];

    X8 = [0.1834346, 0.5255324, 0.7966665, 0.9602899];
    W8 = [0.3626838, 0.3137066, 0.2223810, 0.1012285];

    X9 = [0.0        0.3242534  0.6133714  0.8360311 0.9681602];
    W9 = [0.3302394  0.3123471  0.2606107  0.1806482 0.0812744];

    X10 = [0.1488743 0.4333954 0.6794096 0.8650634 0.9739065];
    W10 = [0.2955242 0.2692667 0.2190864 0.1494513 0.0666713];
    %-----------------------------------------------------------
    % retrieve weights and abscissa
    W = eval(['W' int2str(N)]);   % weight
    X = eval(['X' int2str(N)]);   % abscissa

    r2_r1=r2-r1;
    r3_r1=r3-r1;
    halfarea = rownorm(crossprod(r2_r1,r3_r1))/4;  % size(r1,1) x 1
    evalstr = [F,'(r'];
    for i=1:nargin - 6,
        evalstr = [evalstr,',P',int2str(i)];
    end
    evalstr = [evalstr,')'];

    LEN = size(W,2);
    I = zeros(size(r1,1),Idim);  % Idim is the dimension of integral
    half_r2_r1= r2_r1/2;
    for i=1:LEN,
        t1 = 0.5*(1-X(i));    % scalar
        t2 = 0.5*(1+X(i));
        ct2 = r1+r3_r1*t2;   % common terms
        ct1 = r1+r3_r1*t1;
        s1 =zeros(size(I));
        s2 =zeros(size(I));
        for j=1:LEN,
            d = half_r2_r1*t1;
            r = ct2+ d*(1+X(j));
            Ta = eval(evalstr);
            r = ct2+ d*(1-X(j));
            Tb = eval(evalstr);

            d= half_r2_r1*t2;
            r = ct1+ d*(1+X(j));
            Tc = eval(evalstr);
            r = ct1+ d*(1-X(j));
            Td = eval(evalstr);

            s1 = s1+W(j)*(Ta+Tb);
            s2 = s2+W(j)*(Tc+Td);

        end
        I = I + W(i)*(t1*s1 +t2*s2);
    end

    I = I.*(halfarea*ones(1,Idim));
    %disp('in nit_gq')
    %keyboard

    function nfv = make_surfaces_coincide(fv,Verbose)

    % Align all surfaces in fv to the innermost surface in array patch structure fv.
    % [Align in the sense that resulting vertices are aligned (i.e. coincident)]

    % Centroid of innermost surface
    head_center = mean(fv(1).vertices);
    % Initialize new structure;
    for k=1:length(fv)
        nfv(k) = fv(1);
    end

    for isurf = 2:length(fv) % for each outermost surface
        if Verbose
            makeuswait('start')
            hw = waitbar(0,sprintf('Aligning BEM envelope #%d on innermost envelope #1',isurf));
        end
        for ivert = 1:size(fv(isurf).vertices,1) % project each vertex onto inner surface
            [intersect,indx,t,u,v] = ray_intersect(fv(isurf),fv(1).vertices(ivert,:),-fv(1).vertices(ivert,:)+head_center,'b');
            if isempty(t)
                errordlg(sprintf('Surfaces #%d and #1 are probably not registered into the same coordinate system. Please check registration using BrainStorm''s alignement tool.',isurf)),
                return,
            end
            [tmp, imin] = min(abs(t));
            t = abs(t(imin));

            nfv(isurf).vertices(ivert,:) = intersect(:,imin)';

            if Verbose & ~rem(ivert,50)
                waitbar(ivert/size(fv(isurf).vertices,1),hw)
            end
        end

        %nfv(isurf).faces = fv(1).faces;
        if Verbose
            makeuswait('stop')
            close(hw)
        end

    end

    % -----------------------------------------------------------------

    function fv = make_surfaces_embedded(fv,tol,Verbose)

    % Check that innermost points are within outermost envelopes
    for isurf = 1:length(fv)-1 % for each innermost surface
        iter = 0;
        IS = inpolyhd(fv(isurf).vertices,fv(isurf+1).vertices,fv(isurf+1).faces);
        iWrong = find(IS == 0)'; % index of vertices frominner surface that are either oustide the outersurface (iWrong = 1)
        % or inside but within tol (meters) of outer surface (iWrong  = .5)
        head_center = mean([fv(isurf).vertices]);

        while ~isempty(iWrong) & iter <100 % some supposedly inner points are outside next outermost envelope
            iter = iter+1;   % while on iter avoid infinte loops
            iWrong_out = find(IS == 0)' % Vertices outside surface
            for k = iWrong_out
                [intersect,indx,t,u,v] = ray_intersect(fv(isurf+1),fv(isurf).vertices(k,:),head_center-fv(isurf).vertices(k,:),'b');
                [tmp, imin] = min(abs(t));
                t = t(imin);
                intersect = intersect(:,imin)';
                dist = sign(t)*(fv(isurf).vertices(k,:) - intersect);
                dist = dist/abs(t);
                fv(isurf).vertices(k,:) = intersect - tol*dist;
            end

            IS = inpolyhd(fv(isurf).vertices,fv(isurf+1).vertices,fv(isurf+1).faces);
            iWrong = find(IS==0)'
        end
    end

    % Now test that every inner vertices are at least away of outer surface by tol distance

    % Check that innermost points are within outermost envelopes
    for isurf = 1:length(fv)-1 % for each innermost surface

        head_center = mean([fv(isurf).vertices]);
        if Verbose
            makeuswait('start')
            hw = waitbar(0,sprintf('Checking for minimum separation between surfaces #%d and #%d',isurf, isurf+1));
        end

        for k = 1:size(fv(isurf).vertices,1) % For each vertex
            [intersect,indx,t,u,v] = ray_intersect(fv(isurf+1),fv(isurf).vertices(k,:),head_center-fv(isurf).vertices(k,:),'b');
            % Find closest intersection
            [tmp, imin] = min(abs(t));
            t = t(imin);
            if abs(t) < tol % Vertex is too close: move away from outer surface
                intersect = intersect(:,imin)';
                dist = sign(t)*(fv(isurf).vertices(k,:) - intersect);
                if t == 0
                    t = tol;
                    dist = sign(t)*(fv(isurf).vertices(k,:)+eps - intersect);
                end
                dist = dist/abs(t);
                fv(isurf).vertices(k,:) = intersect - tol*dist;
            end
            if Verbose & ~rem(k,50)
                waitbar(k/size(fv(isurf).vertices,1),hw)
            end
        end
    end
    if Verbose
        makeuswait('stop')
        close(hw)
    end

    % Now that outer surfaces may have moved
    % Check that innermost points of surf#1 are still away enough from surf#2
    for isurf = 1 % for each innermost surface

        head_center = mean([fv(isurf).vertices]);

        for k = 1:size(fv(isurf).vertices,1) % For each vertex
            [intersect,indx,t,u,v] = ray_intersect(fv(isurf+1),fv(isurf).vertices(k,:),head_center-fv(isurf).vertices(k,:),'b');
            % Find closest intersection
            [tmp, imin] = min(abs(t));
            t = t(imin);
            if abs(t) < tol % Vertex is too close: move away from outer surface
                intersect = intersect(:,imin)';
                dist = sign(t)*(fv(isurf).vertices(k,:) - intersect);
                if t == 0
                    t = tol;
                    dist = sign(t)*(fv(isurf).vertices(k,:)+eps - intersect);
                end
                dist = dist/abs(t);
                fv(isurf).vertices(k,:) = intersect - tol*dist;
            end
        end
    end

    % Final check
    % Check that innermost points are within outermost envelopes
    for isurf = 1:length(fv)-1 % for each innermost surface
        iter = 0;
        IS = inpolyhd(fv(isurf).vertices,fv(isurf+1).vertices,fv(isurf+1).faces);
        iWrong = find(IS == 0)' % index of vertices frominner surface that are either oustide the outersurface (iWrong = 1)
        % or inside but within tol (meters) of outer surface (iWrong  = .5)
        head_center = mean([fv(isurf).vertices]);

        while ~isempty(iWrong) & iter <100 % some supposedly inner points are outside next outermost envelope
            iter = iter+1;   % while on iter avoid infinte loops
            iWrong_out = find(IS == 0)' % Vertices outside surface
            for k = iWrong_out
                [intersect,indx,t,u,v] = ray_intersect(fv(isurf+1),fv(isurf).vertices(k,:),-head_center+fv(isurf).vertices(k,:),'o');
                intersect = intersect(:,1)';
                dist = (fv(isurf).vertices(k,:) - intersect);
                dist = dist/abs(t(1));
                fv(isurf).vertices(k,:) = intersect - tol*dist;
            end

            IS = inpolyhd(fv(isurf).vertices,fv(isurf+1).vertices,fv(isurf+1).faces);
            iWrong = find(IS<1)'
        end
    end

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Nov 2 16:17:59 2004
Cross-Directory links are: ON