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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---------------------- 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.
% ------------------------ 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