function [Ginterp] = gain_interp(Rq,bem_gaingrid_mfname)
%GAIN_INTERP - Gain matrix interpolation from a set of pre-computed forward fields defined over grid points in 3D space
% function [Ginterp] = gain_interp(Rq,bem_gaingrid_mfname)
%  FORWARD MODEL CALCULATION USING 3-D INTERPOLATION ON A PRE-COMPUTED GRID
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  FORWARD MODEL CALCULATION USING 3-D INTERPOLATION ON A PRE-COMPUTED GRID
%
% This program computes the (voltage potential)/(radial magnetic field component)
% forward gain matrix for an array of (EEG electrodes)/(MEG magnetometers) 
% interpolating between pre-computed forward model solution points on a densely sampled grid. 
% The grid sampling and associated BEM solution are pre-computed using external programs.
%
% INPUTS: (Required)
%
%  Rq: dipole locations in subjectsc coordinate system (meters)                 (P x 3)
%
%
%  bem_gaingrid_mfname: ".mat" filename which contains the set of 
%                       3-D grid points, variables indicating sampling 
%                       density and the associated pre-computed forward 
%                       model solution.                                    (char string)
%
% OUTPUTS:
%              Ginterp: Forward Gain Solution determined via 3-D 
%                       interpolation between 8 adjacent grid points    (M x 3P)
%
%   John Ermer 03/18/00 (Initial Version)
%   John Ermer 03/21/00 (Optimization of code)
%   John Ermer 03/22/00 (Addition of Persistent Variables)
%   John Ermer 05/21/00 (Added output of Sensor Positions)
%   SB 04/04/2001 (Management of unused entries of the gaingrid matrix )
%   SB 07-03-2003 Major update and code and creation of gain_interp from John Ermer's bem_gain_interp2 
%                 Implementation of inverse-distance weighted interpolation using Shepards method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%<autobegin> ---------------------- 26-May-2004 11:30:18 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Forward Modeling
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\norlig.m
%
% At Check-in: $Author: Mosher $  $Revision: 10 $  $Date: 5/26/04 9:59a $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 24-May-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 26-May-2004 11:30:18 -----------------------



load(bem_gaingrid_mfname);  % Load Pre-computed Fwd Gain Matrix
load(bem_grid_mfname); % load associated grid point locations

GBEM_grid( find(isnan(GBEM_grid(:,1))) ,: ) = [];  % Remove unused entries (for MEG reference channels and EEG reference)

M = size(GBEM_grid,1);             % Number of sensors
P = size(Rq,1);                    % Number of dipole locations

if size(Rq,1) > 100
    hw = waitbar(0,'Forward field interpolation in progress. . .');
end


Ginterp = zeros(size(GBEM_grid,1),3*size(Rq,1));
    
for dip = 1:size(Rq,1)
    if ~rem(dip,200) & exist('hw','var')
        waitbar(dip/size(Rq,1),hw)
    end
    
    % For each dipole - compute trilinear interpolation
    dist2grid = norlig(Rq_bemgrid - repmat(Rq(dip,:), size(Rq_bemgrid,1),1));
    [dist2grid isort] = sort(dist2grid);
    icube = isort(1:8); % Take 8 closest neightbors
    dist2cube = dist2grid(1:8);
    clear dist2grid isort
    %     cube = [Rq_bemgrid(icube,:) - repmat(Rq_bemgrid(icube(1),:),8,1)];
    %     dippos = Rq(dip,:) - Rq_bemgrid(icube(1),:);
    %     dist2cube = norlig(cube - repmat(dippos, 8,1))+eps;

    R = max(dist2cube);
    W = power((R - dist2cube)./ (R*dist2cube),2); % Interpolation weights from Shepards method
    W = W/sum(W);

    for k = 1:length(W)-1 % length(W)-1 because W(end) = 0 by construction.
        Ginterp(:,3*(dip-1)+1) = Ginterp(:,3*(dip-1)+1) + W(k) * GBEM_grid(:,3*(icube(k)-1)+1) ;
        Ginterp(:,3*(dip-1)+2) = Ginterp(:,3*(dip-1)+2) + W(k) * GBEM_grid(:,3*(icube(k)-1)+2) ;
        Ginterp(:,3*(dip-1)+3) = Ginterp(:,3*(dip-1)+3) + W(k) * GBEM_grid(:,3*(icube(k)-1)+3) ;
    end
    
end

if exist('hw','var')
    close(hw)
end

