[Master Index] [Index for Toolbox]

# bem_kernel

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

### Function Synopsis

    [Kv,Kb]= bem_kernel(rQ,basis,test,geometry,nodes,cdv,Te,P_wts,Tm,S);

### Help Text

    BEM_KERNEL - computes kernel matrices for EEG and MEG
     function [Kv,Kb]= bem_kernel(rQ,basis,test,geometry,nodes,cdv,Te,P_wts,Tm,S);
    *******************************************************************************

     This computes kernel matrices for EEG and MEG
           Gv=Te*Gq
           Gb= Gp+Tm*Gq
           Te and Tm are computed off-line using bem1.m

      Input parameters:
      rQ        : dipole locations,                            no_dipole x 3
      To get the following, just load the .mat file saved by bem1.m
      basis     : string --  constant bem or linear bem
      test      : string --  collocation, galerkin
      N         : [no of EEG sensors, total no of triangles/ nodes]   2 x1
                  e.g. [492,1476].
      geometry,nodes,cdv : see bem1.m for description
      S         : MEG sensor locations
                       no_sensor x 3
      P_wts     : EEG Triangle Node Weights (Applies to Linear Basis only;
                  (Set to empty if not used)                          (Meeg x 3)
      Te,Tm     : transfer matrices
      Be sure to provide ALL input parameters

      Output Parameters:
      Kv        : kernel matrix of EEG ,              no_eegsensor x (3*no_dipoles)
                  Kv(i,3*j-2:3*j): EEG due to the x,y,and z component of the j-th
                  dipole at the i-th EEG sensor

      Kb        : kernel matrix of MEG ,              (3*no_megsensor) x (3*no_dipole)
                  Kb(3*i-2:3*i,3*j-2):  MEG due to x component of the j-th dipole
                  Kb(3*i-2:3*i,3*j-1):  MEG due to y component of the j-th dipole
                  Kb(3*i-2:3*i,3*j)  :  MEG due to z component of the j-th dipole

      Modifications:
                 - 040200: Optimized code for Constant Galerkin Option (J. Ermer)
                 - 042800: Added Logic to interpolate between triangle nodes for
                           linear basis (J. Ermer)
                 - SB 07-Mar-2003 : moved all Ermer's external functions as subfunctions

    ******************************************************************************
    profile on -detail builtin

    %%% This part checks validity of input parameters %%%%

### Cross-Reference Information

This function calls

    crossprod                          C:\BrainStorm_2001\Toolbox\crossprod.m
    dotprod                            C:\BrainStorm_2001\Toolbox\dotprod.m
    norlig                             C:\BrainStorm_2001\Toolbox\norlig.m
    rownorm                            C:\BrainStorm_2001\Toolbox\rownorm.m

This function is called by

    bem_gain                           C:\BrainStorm_2001\Toolbox\bem_gain.m

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

    function [Kv,Kb]= bem_kernel(rQ,basis,test,geometry,nodes,cdv,Te,P_wts,Tm,S);
    %BEM_KERNEL - computes kernel matrices for EEG and MEG
    % function [Kv,Kb]= bem_kernel(rQ,basis,test,geometry,nodes,cdv,Te,P_wts,Tm,S);
    %*******************************************************************************
    %
    % This computes kernel matrices for EEG and MEG
    %       Gv=Te*Gq
    %       Gb= Gp+Tm*Gq
    %       Te and Tm are computed off-line using bem1.m
    %
    %  Input parameters:
    %  rQ        : dipole locations,                            no_dipole x 3
    %  To get the following, just load the .mat file saved by bem1.m
    %  basis     : string --  constant bem or linear bem
    %  test      : string --  collocation, galerkin
    %  N         : [no of EEG sensors, total no of triangles/ nodes]   2 x1
    %              e.g. [492,1476].
    %  geometry,nodes,cdv : see bem1.m for description
    %  S         : MEG sensor locations
    %                   no_sensor x 3
    %  P_wts     : EEG Triangle Node Weights (Applies to Linear Basis only;
    %              (Set to empty if not used)                          (Meeg x 3)
    %  Te,Tm     : transfer matrices
    %  Be sure to provide ALL input parameters
    %
    %  Output Parameters:
    %  Kv        : kernel matrix of EEG ,              no_eegsensor x (3*no_dipoles)
    %              Kv(i,3*j-2:3*j): EEG due to the x,y,and z component of the j-th
    %              dipole at the i-th EEG sensor
    %
    %  Kb        : kernel matrix of MEG ,              (3*no_megsensor) x (3*no_dipole)
    %              Kb(3*i-2:3*i,3*j-2):  MEG due to x component of the j-th dipole
    %              Kb(3*i-2:3*i,3*j-1):  MEG due to y component of the j-th dipole
    %              Kb(3*i-2:3*i,3*j)  :  MEG due to z component of the j-th dipole
    %
    %
    %  Modifications:
    %             - 040200: Optimized code for Constant Galerkin Option (J. Ermer)
    %             - 042800: Added Logic to interpolate between triangle nodes for
    %                       linear basis (J. Ermer)
    %             - SB 07-Mar-2003 : moved all Ermer's external functions as subfunctions
    %
    %******************************************************************************
    %profile on -detail builtin
    %
    %%%% This part checks validity of input parameters %%%%
    %

    %<autobegin> ---------------------- 14-Jun-2004 17:09:46 -----------------------
    % --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
    %
    % CATEGORY: Forward Modeling
    %
    % Alphabetical list of external functions (non-Matlab):
    %   toolbox\crossprod.m
    %   toolbox\dotprod.m
    %   toolbox\norlig.m
    %   toolbox\rownorm.m
    %
    % Subfunctions in this file, in order of occurrence in file:
    %   I = nit_gq_lg(tri_data,Rq)
    %   g = glg1(r,rq,rnxrn1)
    %   I = nit_gq_cg(r1,r2,r3,Rq)
    %   g = gterm_constant(r,rq)
    %
    % At Check-in: $Author: Mosher $  $Revision: 15 $  $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:46 -----------------------

    if size(rQ,2)~=3,
        error(' rQ must have three columns')
    end
    %
    %%%% This part determines BEM type and EEG and/or MEG Processing %%%%
    %
    if ~isempty(Te), mode =1; end
    if ~isempty(Tm), mode=2; end
    if ~isempty(Te) & ~isempty(Tm), mode=3; end

    if lower(basis(1)) =='c', basis_opt = 0; end  % constant BEM
    if lower(basis(1)) =='l', basis_opt =1; end  % linear BEM

    if lower(test(1)) =='c', choice =0; end  % collocation
    if lower(test(1)) =='g', choice =1; end  % galerkin
    %
    %%%% This part preallocates variables and computes triangle vertices %%%%
    %
    no_dipole =  size(rQ,1);
    %
    r1=nodes(geometry(:,1),1:3);  % vertices of the ith triangle
    r2=nodes(geometry(:,2),1:3);
    r3=nodes(geometry(:,3),1:3);
    %
    %Gq=zeros(N(2),3*no_dipole);
    if ~isempty(Te)
        Gq=zeros(size(Te,2),3*no_dipole);
    else
        Gq=zeros(size(Tm,2),3*no_dipole);
    end

    %
    %%% This part performs processing unique to Constant and Linear Basis Functions %%%%
    %
    if basis_opt == 0,                          % Constant BEM
        whichsurf = nodes(geometry(:,1),4);
        r = (r1+r2+r3)/3;                         % centroid of each triangle
        no_elt = size(geometry,1);
    else                                        % Linear BEM
        whichsurf = nodes(:,4);
        r = nodes(:,1:3);                         % position of each node
        no_node = size(nodes,1);
    end
    %
    cdvsum = cdv(whichsurf,1)+cdv(whichsurf,2);    % conductivity sum (size depends
    %                                                on Basis Function used
    %%%% This part computes Gq for Each of the BEM Methods %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    if choice ==0,                              % collocation (constant -or- linear)
        c = 2*pi*cdvsum;                         % (vectorized)
        Gq = gterm_constant(r,rQ);
        Gq = Gq./repmat(c,1,3*no_dipole);
        %
    else                                        % Galerkin forms
        %
        if basis_opt ==0,                         % Constant Galerkin (vectorized)
            %
            %      c = 2*pi*cdvsum.*area;
            c = 4*pi*cdvsum;        % (Residual term after area/2 cancels out in nit_Gq_cg)
            Gq = nit_gq_cg(r1,r2,r3,rQ);
            Gq = Gq./repmat(c,1,3*no_dipole);
            %
        elseif basis_opt == 1,                   % Linear Galerkin, 4x4 point GQ
            %
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            whichsurf = nodes(geometry(:,1),4);
            cdvsum = cdv(whichsurf,1)+cdv(whichsurf,2);    % conductivity sum
            %
            DET = dotprod(r1,crossprod(r2,r3));  % Nt x 1
            rnxrn1 = [crossprod(r2,r3),crossprod(r3,r1),crossprod(r1,r2)]; % (Nt x 9)
            %
            ptri = zeros(no_node,3);
            tri_data = zeros(6*no_node,15);
            t_cnt = 0;
            %
            for n=1:no_node  % Loop through all nodes and develop "triangle list"
                ptri = [];
                for  i=1:3,    % find triangles containing this node
                    t =find(geometry(:,i)==n);
                    ptri= [ptri; t,i*ones(size(t))];  % Save record of all triangles
                end
                %
                for i=1:size(ptri,1)  % Loop through all triangles attached to node
                    t= ptri(i,1);     % triangle number
                    t_cnt = t_cnt + 1; % Running count of all triangles
                    tri_data(t_cnt,:) =[n t ptri(i,2) r1(t,:) r2(t,:) r3(t,:) ....
                            rnxrn1(t,3*ptri(i,2)-2:3*ptri(i,2))];
                end
            end  % node index
            tri_data = tri_data(1:t_cnt,:);  % Delete excess terms
            %
            %%% This part integrates and performs post-scaling on all triangles %%%%
            %
            Gq_all = nit_gq_lg(tri_data,rQ);  % Integrate all triangles
            term = (2*pi)*cdvsum(tri_data(:,2),1).*DET(tri_data(:,2),1);
            Gq_all = Gq_all./repmat(term,1,3*no_dipole);
            %
            %%%% This part sums together all triangles associated with a given node %%%%
            %
            for n=1:no_node
                indx = find(tri_data(:,1)==n);
                Gq(n,:) = sum(Gq_all(indx,:));
            end
            %
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        end
    end
    %
    %%%% This part computes the final Kernel for EEG %%%%
    %
    Kv=[];
    if mode==1 | mode==3,
        if basis_opt==1  % Linear Basis (Interpolated Between 3 Triangle Vertices)
            Kv = zeros(size(Te,1)/3,3*no_dipole);
            for m=1:(size(Te,1)/3)
                Kv(m,:) = sum(repmat(P_wts(m,:)',1,3*no_dipole).*(Te(3*m-2:3*m,:)*Gq));
            end
        else
            Kv = Te*Gq;   % Constant Basis
        end
    end
    %
    %%%% This part computes the final Kernel for MEG %%%%
    %
    Kb=[];
    if mode==2 | mode ==3
        no_sensor=size(S,1);
        m=3*no_sensor;
        Gp=zeros(m,3*no_dipole);
        % compute Gp
        for i=1:no_dipole
            S_rQ=[S(:,1)-rQ(i,1),S(:,2)-rQ(i,2),S(:,3)-rQ(i,3)];
            mag3=rownorm(S_rQ).^3;
            t= S_rQ ./[mag3,mag3,mag3];
            Gp(1:3:m,[3*i-1,3*i])=[t(:,3),-t(:,2)];
            Gp(2:3:m,[3*i-2,3*i])=[-t(:,3),t(:,1)];
            Gp(3:3:m,[3*i-2,3*i-1])=[t(:,2),-t(:,1)];
        end

        Kb=Gp+Tm*Gq;      % u0/4pi OMITTED HERE
    end

    % SUBFUNCTIONS -----------------------------------------------------------------------

    function I = nit_gq_lg(tri_data,Rq)
    %NIT_GQ_LG : numerical integration over tessel tri's for const galerkin function
    % function I = nit_gq_lg(tri_data,Rq)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % This function integrates the vector function gterm_constant
    % (for the constant galerkin option function over triangles
    % with vertices r1,r2,and r3 using 4-point Gauss Quadrature
    %
    % where  r=(x,y,z) is a point on the triangle, and P1,P2,.... are parameters
    % defined in F.
    %
    % Inputs:
    % r1, r2,r3 : specify the vertices of triangles over which function is to be
    %             integrated (Nt x 3)
    %       Rq  : Dipole Positions (P x 3)
    %
    % Output:
    %         I : integral (Nt x 3P)
    %
    % %%% John Ermer 04/02/00
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %<autobegin> -------- 20-Nov-2002 14:07:25 ------------------------------
    % ---- Automatically Generated Comments Block using auto_comments -----------
    %
    % Alphabetical list of external functions (non-Matlab):
    %   toolbox\crossprod.m
    %   toolbox\glg1.m
    %   toolbox\rownorm.m
    %
    % At Check-in: $Author: Mosher $  $Revision: 15 $  $Date: 6/14/04 3:37p $
    %
    % Copyright (c) 2002 BrainStorm MMII
    % Source code may not be distributed in original or altered form.
    % See bst_splashscreen, http://neuroimage.usc.edu, or email leahy@sipi.usc.edu
    %   for license and copyright notices.
    %<autoend> ---------- 20-Nov-2002 14:07:25 ------------------------------

    %profile on -detail builtin
    %
    r1 = tri_data(:,4:6);
    r2 = tri_data(:,7:9);
    r3 = tri_data(:,10:12);
    rnxrn1 = tri_data(:,13:15);
    %
    %-----------------------------------------------------------
    % W: weight table, X : abscissa table
    X = [0.3399810, 0.8611363];
    W = [0.6521452, 0.3478548];
    %-----------------------------------------------------------
    %
    Nt = size(r1,1);
    P = size(Rq,1);
    r2_r1=r2-r1;   % (Nt x 3)
    r3_r1=r3-r1;   % (Nt x 3)
    halfarea = rownorm(crossprod(r2_r1,r3_r1))/4;  % (Nt x 1)
    %
    I = zeros(Nt,3*P);    % dimension of integral is 3
    half_r2_r1= r2_r1/2;  % (Nt x 3)
    %
    %%% This part loops through 16 points on triangle, evaluates them, and sums them
    %
    for i=1:2
        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:2
            d = half_r2_r1*t1;
            r = ct2+ d*(1+X(j));
            Ta = glg1(r,Rq,rnxrn1);
            r = ct2+ d*(1-X(j));
            Tb = glg1(r,Rq,rnxrn1);
            %
            d= half_r2_r1*t2;
            r = ct1+ d*(1+X(j));
            Tc = glg1(r,Rq,rnxrn1);
            r = ct1+ d*(1-X(j));
            Td = glg1(r,Rq,rnxrn1);
            %
            s1 = s1+W(j)*(Ta+Tb);
            s2 = s2+W(j)*(Tc+Td);
            %
        end
        I = I + W(i)*(t1*s1 +t2*s2);
    end
    %
    I = I.*repmat(halfarea,1,3*P);

    % Subfunctions%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function g = glg1(r,rq,rnxrn1)
    %GLG1 (overwrite succinct one line summary here)
    % function g = glg1(r,rq,rnxrn1)
    % gterm_constant: Integration Function for Linear Galerkin Method
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % Inputs:
    %         r : Location Points on Triangles (Nt x 3)
    %        rq : Dipole Positions (P x 3)
    %     rnxrn1: cross-product term (Ntx3
    %
    % Output:
    %         g : Function Result (Nt x 3P)
    %
    % %%% John Ermer 04/02/00
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %<autobegin> -------- 20-Nov-2002 14:05:57 ------------------------------
    % ---- Automatically Generated Comments Block using auto_comments -----------
    %
    % Alphabetical list of external functions (non-Matlab):
    %   toolbox\dotprod.m
    %   toolbox\rownorm.m
    %
    % At Check-in: $Author: Mosher $  $Revision: 15 $  $Date: 6/14/04 3:37p $
    %
    % Copyright (c) 2002 BrainStorm MMII
    % Source code may not be distributed in original or altered form.
    % See bst_splashscreen, http://neuroimage.usc.edu, or email leahy@sipi.usc.edu
    %   for license and copyright notices.
    %<autoend> ---------- 20-Nov-2002 14:05:57 ------------------------------

    %
    Nt = size(r,1);
    P = size(rq,1);
    %
    %%%% Modified version incorporating vectorized processing
    %
    r_rq = repmat(r,P,1) - reshape(repmat(rq',Nt,1),3,Nt*P)'; % (Nt*P)x(3)
    n3 = norlig(r_rq)'; %rownorm(r_rq);                                       % (Nt*P)x(3)
    n3 = n3.*n3.*n3;                                          % (Nt*P)x(3)
    dprod = repmat(dotprod(rnxrn1,r),1,P);
    %
    g = zeros(Nt,3*P);
    g(:,[1:3:3*P-2]) = dprod.*reshape(r_rq(:,1)./n3,Nt,P);  % x-comp
    g(:,[2:3:3*P-1]) = dprod.*reshape(r_rq(:,2)./n3,Nt,P);  % y-comp
    g(:,[3:3:3*P]) = dprod.*reshape(r_rq(:,3)./n3,Nt,P);    % z-comp
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %

    function I = nit_gq_cg(r1,r2,r3,Rq)
    %NIT_GQ_CG : numerical integration over tessel tri's for const galerkin function
    % function I = nit_gq_cg(r1,r2,r3,Rq)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % This function integrates the vector function gterm_constant
    % (for the constant galerkin option function over triangles
    % with vertices r1,r2,and r3 using 4-point Gauss Quadrature
    %
    % where  r=(x,y,z) is a point on the triangle, and P1,P2,.... are parameters
    % defined in F.
    %
    % Inputs:
    % r1, r2,r3 : specify the vertices of triangles over which function is to be
    %             integrated (Nt x 3)
    %       Rq  : Dipole Positions (P x 3)
    %
    % Output:
    %         I : integral (Nt x 3P)
    %
    % %%% John Ermer 04/02/00
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %<autobegin> -------- 20-Nov-2002 14:07:23 ------------------------------
    % ---- Automatically Generated Comments Block using auto_comments -----------
    %
    % Alphabetical list of external functions (non-Matlab):
    %   toolbox\gterm_constant.m
    %
    % At Check-in: $Author: Mosher $  $Revision: 15 $  $Date: 6/14/04 3:37p $
    %
    % Copyright (c) 2002 BrainStorm MMII
    % Source code may not be distributed in original or altered form.
    % See bst_splashscreen, http://neuroimage.usc.edu, or email leahy@sipi.usc.edu
    %   for license and copyright notices.
    %<autoend> ---------- 20-Nov-2002 14:07:23 ------------------------------

    %profile on -detail builtin
    %
    %-----------------------------------------------------------
    % W: weight table, X : abscissa table
    X = [0.3399810, 0.8611363];
    W = [0.6521452, 0.3478548];
    %-----------------------------------------------------------
    %
    Nt = size(r1,1);
    P = size(Rq,1);
    r2_r1=r2-r1;   % (Nt x 3)
    r3_r1=r3-r1;   % (Nt x 3)
    %
    I = zeros(Nt,3*P);    % dimension of integral is 3
    half_r2_r1= r2_r1/2;  % (Nt x 3)
    %
    %%% This part loops through 16 points on triangle, evaluates them, and sums them
    %
    for i=1:2
        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:2
            d = half_r2_r1*t1;
            r = ct2+ d*(1+X(j));
            Ta = gterm_constant(r,Rq);
            r = ct2+ d*(1-X(j));
            Tb = gterm_constant(r,Rq);
            %
            d= half_r2_r1*t2;
            r = ct1+ d*(1+X(j));
            Tc = gterm_constant(r,Rq);
            r = ct1+ d*(1-X(j));
            Td = gterm_constant(r,Rq);
            %
            s1 = s1+W(j)*(Ta+Tb);
            s2 = s2+W(j)*(Tc+Td);
            %
        end
        I = I + W(i)*(t1*s1 +t2*s2);
    end
    %
    %I = I.*repmat(halfarea,1,3*P); % Cancels out in calling program

    function g = gterm_constant(r,rq)
    %gterm_constant
    % function g = gterm_constant(r,rq)

    %<autobegin> -------- 20-Nov-2002 14:06:02 ------------------------------
    % ---- Automatically Generated Comments Block using auto_comments -----------
    %
    % Alphabetical list of external functions (non-Matlab):
    %   toolbox\rownorm.m
    %
    % At Check-in: $Author: Mosher $  $Revision: 15 $  $Date: 6/14/04 3:37p $
    %
    % Copyright (c) 2002 BrainStorm MMII
    % Source code may not be distributed in original or altered form.
    % See bst_splashscreen, http://neuroimage.usc.edu, or email leahy@sipi.usc.edu
    %   for license and copyright notices.
    %<autoend> ---------- 20-Nov-2002 14:06:02 ------------------------------

    if size(rq,1)  == 1 % Just one dipole
        r_rq= [r(:,1)-rq(1),r(:,2)-rq(2),r(:,3)-rq(3)];
        n = rownorm(r_rq).^3;
        g = r_rq./[n,n,n];
    else
        g = zeros(size(r,1),3*size(rq,1));
        isrc = 1;
        for k = 1:size(rq,1)
            r_rq= [r(:,1)-rq(k,1),r(:,2)-rq(k,2),r(:,3)-rq(k,3)];
            n = rownorm(r_rq).^3;
            g(:,3*(isrc-1)+1: 3*isrc) = r_rq./[n,n,n];
            isrc = isrc + 1;
        end

    end

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