function HeadModel = source_grids(HeadModel,Channel,SubjectFV,GUI);
%SOURCE_GRIDS - Generate source grids for parametric model searches
%function HeadModel = source_grids(HeadModel,Channel,SubjectFV,GUI);
%
% HeadModel and Channel are the standard structures:
% HeadModel:
%  .Function
%  .Param, struct, e.g. 'center','radius'
%  .GridLoc, RETURNED: cell array generated by this program
%  .SearchTess, unused, maybe dropped
%  .Gain, RETURNED: cell array one-to-one with GridLoc
%  .GridName
%  .GridLoc
%  .ImageGridOrient
%  .Gain
%  .GainCovar
%
% Because of the size of the Image Grids, I suggest that .Image* fields NOT 
%  be sent along in the calling structure, but added in later.
%
% Channel:
%  .Loc
%  .Orient
%  .Weight
%  .Type
%  .Comment
%
% SubjectFV is the faces and vertices from the SubjectTess structure. Note the 
%  structures specific to Matlab, not BrainStorm: lowercase names and row orientation.
% SubjectFV.faces is n x 3, used only to reducepatch the vertices to around 250.
% SubjectFV.vertices is 3 x m, used to bound the search grid, should 
%  be a smooth inner-skull-like surface inside of which the points will be generated.
%
% GUI : a structure passing several parameters to the source_grids
% .VALIDORDER = Source orders Handled by the forward models.
%       As of 03/13/00 : [-1 0 1] in MEG and [-1] for EEG
%       Initially non-synchronous sources only
% .SPACING = ; initial grid spacing in mm for the SearchGrid 
% .VERBOSE : 1 (res. 0) turns on (res. off) the verbose mode
% .MEG and .EEG = two binary flags indicating whether the grid gain needs to be computed (1) for either modality or not (0) 
%

%<autobegin> ---------------------- 02-Nov-2004 16:16:30 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Forward Modeling
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bst_message_window.m
%   toolbox\closest.m
%   toolbox\colnorm.m
%   toolbox\dip.m
%   toolbox\good_channel.m
%   toolbox\rownorm.m
%   toolbox\subcorr_gui.m
%
% At Check-in: $Author: Mosher $  $Revision: 21 $  $Date: 11/02/04 2:38p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 02-Nov-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> ------------------------ 02-Nov-2004 16:16:30 -----------------------


% /---Script Authors--------------------------------------\
% |                                                      |
% | *** John C. Mosher, Ph.D.                            |
% |  Biophysics Group                                    |
% |                                                      |
% | *** Sylvain Baillet Ph.D.                            |
% | Cognitive Neuroscience & Brain Imaging Laboratory    |
% | CNRS UPR640 - LENA                                   | 
% | Hopital de la Salpetriere, Paris, France             |
% | sylvain.baillet@chups.jussieu.fr                     |
% |                                                      |
% \------------------------------------------------------/


% Script History -----------------------------------------------------------------------------------------------------------
% Date of creation: October, 25 1999
%
% Nov 2001: Checked sequential computation of EEG and MEG head models when requested
% SB March 2002 :  - Added VERBOSE mode 
%                  - Can handle simulateous creation of MEG and EEG GridLoc and SearchGridGain           
% SB 24-Apr-2002 : - Script header platic surgery
%                  - Moved back to .Faces and .Vertices field names for SubjectFV
% SB 26-Apr-2002 : - Moved to BsT MMII new HeadModel filednames (e.g. .GridName)
%                  - Removed computation of Gains for CD or CME pairs : will sample them from original CD and CME Gains (like we already did actually)
%                  - Added the GUI.MEG and GUI.EEG flags
% SB 01-Jul-2002 : - Added error detection when EEG channels are used for EOG or ECoG.
% SB 02-Jul-2002 : - Fixed bug when gridding for CME only
% SB 12-Sep-2002 : - SubjectFV is a tessellation structure that depicts the inner-skull volume
%                    where the search volumic grid must be contained
%                    Former version used the scalp surface which triggered errors in the 3-shell Legendre gain matrix computation
%                    because it checked for sources outside the inner sphere. 
% SB 19-Nov-2002 : - Minor display upgrade of display
%                    Updated management of EEG reference.
% JCM 20-Nov-2002    Minor fix of duplicate autocomments blocks, some comment formatting
% JCM 18-Nov-2003    Fixing problem of vertice being m x 3 or 3 x m. Force to be 3 x m.
% JCM 19-Oct-2004    Increased CORR hardwire to be 0.98 for tighter spacing
% ------------------------------------------------------------------


SPACING = GUI.SPACING/1000; % in meters
VALIDORDER = GUI.VALIDORDER;
SourceOrderString = {'Current Dipole',...
      'Unvalid Order (was Magnetic Dipole Moment)',...
      'Current Multipole Expansion',...
      'Current Dipole Pairs',...
      'Unvalid Order (was Magnetic Dipole Moment PAIRS)',...    
      'Current Multipole Expansion Pairs'};

% Local hardwired options
SENSOR_PROXIMITY = 25/1000; % grid can't be closer than this to a sensor coil

% pick a pruning technique 
%PRUNE = 'Power'; % based on the remnant power in the projected gain matrices
PRUNE = 'Subspace'; % based on the correlation of adjacent points.
VERBOSE = 0;
% Pruning Options
LOCAL = SPACING*2; % local sphere to check pruning options
CORR = 0.98; % grid points must be less correlated than this to remain good


[m,n] = size(SubjectFV.vertices);

if m ~= 3,
   SubjectFV.vertices = SubjectFV.vertices'; % should be 3 x n, each column a vector
end

% MEG & EEG respective channel indices
% Find MEG and EEG CHANndx
ChannelFlag = ones(size(Channel)); 
MEGndx = good_channel(Channel,ChannelFlag,'MEG');
EEGndx = good_channel(Channel,ChannelFlag,'EEG');


HeadModel.GridLoc = cell(1,length(VALIDORDER));
HeadModel.Gain = cell(1,length(VALIDORDER));

% Check valid source orders depeding on the requests
if ~GUI.MEG | (GUI.MEG & strcmp(HeadModel.Function{MEGndx(1)},'meg_bem')) % 'doing MEG BEM or EEG only -> Authorize current dipole models
   inotvalid = find(...
      VALIDORDER == 0 | ...
      VALIDORDER == 1 | ...
      VALIDORDER == 3 | ...
      VALIDORDER == 4);
   
   VALIDORDER(inotvalid) = []; % remove all unvalid source configurations
   
end

for ORDER = VALIDORDER, %the source model orders
   iOrder = ORDER + 2; % associated indexer, starting with 1.
   
   if 0,%GUI.VERBOSE
      bst_message_window({sprintf('Creating %s gain matrices for volume source scan ', SourceOrderString{iOrder})}),
   end
   
   switch ORDER % now translate order to model orders
   case {-1,0,1}
      % order remains the same
      SYNC = 1; % single source
   case {2,3,4}
      ORDER = ORDER - 3; % revert to true order
      SYNC = 2; % dual synchronous source
   end
   
   % We anticipate that single dipoles need the densest grid. Each subsequent
   %  order uses the previous grid.
   
   if ORDER==VALIDORDER(1), % first time through
      % first generate a plaid grid bounded by the subject vertices
      
      % only need about 250 vertices to box in the grid.
      if size(SubjectFV.faces,1)>500
         tempFV = SubjectFV;
         tempFV.vertices = tempFV.vertices'; % so that they are n x 3
         NFV = reducepatch(tempFV,500); % number of faces
         GridBox = NFV.vertices'; % 3 x ~500
      else
         GridBox = SubjectFV.vertices;
      end
      
      mx = max(GridBox,[],2);
      mn = min(GridBox,[],2);
      
      [X,Y,Z] = meshgrid([mn(1):SPACING:mx(1)],...
         [mn(2):SPACING:mx(2)],[mn(3):SPACING:mx(3)]); 
      
      L = [X(:) Y(:) Z(:)]';
      tempLn = rownorm(L);
      L(find(~tempLn),:) = []; % delete points at the origin as degenerate in MEG
      
      
      % prune sources outside the the GridBox
      % We find the vertex closest to the grid point. If the radius of that
      %  grid point is greater than the vertex, we mark it bad.
      % Inelegant 3-d measure of "inside", could use better.
      % "radius" is defined to a central mean of the vertices.
      mnVert = mean(SubjectFV.vertices,2); % mean of the original vertices.
      
      Good = ones(1,size(L,2));   
      if GUI.VERBOSE
         hwait = waitbar(0,sprintf('Sifting %.0f dipoles plaid at %.1f mm',...
            size(L,2),SPACING*1000));
      end
      
      for i = 1:size(L,2),
         if(~rem(i,100)),
            if GUI.VERBOSE, waitbar(i/size(L,2),hwait); end
         end
         pt2 = closest(L(:,i)',GridBox'); % old routine using rows of data not cols.
         if(norm(L(:,i)-mnVert) > 0.99*norm(pt2'-mnVert)),
            Good(i) = 0; % mark it bad
         end
      end
      if GUI.VERBOSE
         close(hwait), end
      
      if GUI.MEG & GUI.EEG
         ChanNdx = MEGndx;
      elseif GUI.MEG
         ChanNdx = MEGndx;
      elseif GUI.EEG 
         ChanNdx = EEGndx;            
      end
      
      % should also prune out stuff adjacent sensors.
      if GUI.VERBOSE
         hwait = waitbar(0,sprintf('Checking %.0f Sensors for %.1f mm Proximity',...
            length(Channel(ChanNdx)),SENSOR_PROXIMITY*1000));
      end
      
      for i =1:length(Channel(ChanNdx)),
         if(~rem(i,50)),
            if GUI.VERBOSE, waitbar(i/length(Channel(ChanNdx)),hwait); end
         end      
         for j = 1:size(Channel(ChanNdx(1)).Loc,2), % each coil
            DIST = colnorm(Channel(ChanNdx(i)).Loc(:,j+zeros(1,size(L,2))) - L);
            ndx= find(DIST < SENSOR_PROXIMITY); % don't get closer than this
            Good(ndx) = 0; % mark bad
         end
      end
      if GUI.VERBOSE, close(hwait), end
      Gndx = find(Good); % the good grid points
      L = L(:,Gndx); % downselect
   else
      % next order, use the previous grid to start
      L = HeadModel.GridLoc{iOrder-2}; % previous grid (Synchronous sources)    
   end
   
   G = cell(1,2);
   
   if xor(GUI.MEG,GUI.EEG)
      
      if GUI.MEG
         ChanNdx = MEGndx;
         sstring = 'MEG';
      else
         ChanNdx = EEGndx;
         sstring = 'EEG';
      end
      
      if GUI.VERBOSE
         bst_message_window(sprintf('Creating %s gain matrix for %.0f sensors and %.0f sources...',...
            sstring, length(ChanNdx),size(L,2))); 
      end
      
      %G = feval(HeadModel.Function{ChanNdx(1)},L,Channel(ChanNdx),HeadModel.Param(ChanNdx),ORDER); % same for synchronous sources as well
      G = feval(HeadModel.Function{ChanNdx(1)},L,Channel,HeadModel.Param,ORDER); % same for synchronous sources as well
      
   else % Both MEG and EEG are requested 
      if GUI.VERBOSE
         bst_message_window(sprintf('Creating MEG gain matrix first, for %.0f sensors and %.0f sources...',...
            length(MEGndx),size(L,2))); 
      end
      
      %G = feval(HeadModel.Function{MEGndx(1)},L,Channel(MEGndx),HeadModel.Param(MEGndx),ORDER); % same for synchronous sources as well    
      G = feval(HeadModel.Function{MEGndx(1)},L,Channel,HeadModel.Param,ORDER); % same for synchronous sources as well    
      
   end
   
   DIM = size(G,2)/size(L,2);% dimension of each source
   
   % done generating initial grid and gain matrix
   
   % now prune the grid of redundant source points
   
   Good = ones(1,size(L,2)); % every grid point initialized as good
   Goodndx = reshape([1:size(G,2)]',DIM,size(G,2)/DIM); % indexer into gain matrix, by columns
   % therefore Goodndx(DIM,:)/DIM is an indexer into Good.
   
   PpG = zeros(size(G)); % zero out all locations
   if(strcmp(PRUNE,'Power')),
      % need a knee calculation, where NOS is the std dev of the white noise
      KNEE = NOS*sqrt(size(G,1)); % sum of the total sensor noise
   end
   if GUI.VERBOSE
      hwait = waitbar(0,sprintf('Processing %.0f sources using %s',...
         length(Good),PRUNE)); end
   
   for i = 1:length(Good),
      NumGoodSources = sum(Good); % how many good sources are there still left
      if(Good(i)), % if it's still good, i.e. hasn't been deleted
         
         dndx = [-(DIM-1):0]+i*DIM; % indexer to gain matrix
         
         % Check only those sources in geometric proximity to the source.
         %  Done for speed considerations.
         DIST = colnorm(L(:,i+zeros(1,NumGoodSources)) - L(:,find(Good)));
         localndx = find(DIST <= LOCAL);
         if 0%GUI.VERBOSE
            bst_message_window(sprintf('Processing %.0f with %.0f local out of %.0f good sources. . .',...
               i,length(localndx),NumGoodSources));            end
         
         switch PRUNE
         case 'Power'
            
            Ug = orth(G(:,dndx)); % orthogonalize the source
            % project source matrices away from the present point.
            % but only bother with the good points. The other columns will be left alone.
            %PpG = G - Ug*(Ug'*G);
            PpG(:,Goodndx(:,localndx)) = G(:,Goodndx(:,localndx)) - ...
               Ug*(Ug'*G(:,Goodndx(:,localndx)));
         case 'Subspace'
            % do nothing
         end
         
         if GUI.VERBOSE & (~rem(i,100))
            waitbar(i/length(Good),hwait); end
         
         for j = localndx(end:-1:1), % for all remaining good sources, in reverse search order
            Source_Num = Goodndx(DIM,j)/DIM; % the source number
            switch PRUNE
            case 'Power'
               Strongest_Source = svd(PpG(:,Goodndx(:,j))); % decompose the projected matrix
               Strongest_Source = Strongest_Source(1); % the strongest
               Strongest_Source = Strongest_Source * DIP; % scaled by anticipated dipole strength
               
               % now standard deviation of the noise is NOS, the "knee" in the equation.
               % if the Strongest_Source is below the knee, delete it.
               
               if(Strongest_Source < KNEE), % weaker than noise
                  if(Source_Num ~= i), % but don't delete the source itself
                     if 0%GUI.VERBOSE
                        bst_message_window(sprintf('Deleting Source %.0f',Source_Num));
                     end          
                     Good(Source_Num) = 0; % mark it deleted
                     Goodndx(:,j) = []; % remove this column
                     % because of reverse sequencing, has no effect on the lower indices in the for loop
                  end
               end
            case 'Subspace'
               nextCorr= subcorr_gui(G(:,dndx),G(:,Goodndx(:,j)),...
                  struct('REG','Tikhonov','Tikhonov',100));
               if isempty(nextCorr) % Occurs when G is full of zeros. Happens for electrodes wrongly involved in forward modeling though they were measuring non-EEG data (like EOG or ECoG)
                  errordlg('Please make sure sensor coordinates are all correct and non-zero','Forward modeling aborted')
                  return
               end
               
               if(nextCorr(end) > CORR), % all corrs too close in correlations
                  if(Source_Num ~= i), % but don't delete the source itself
                     if 0%GUI.VERBOSE,
                        bst_message_window(sprintf('Deleting Source %.0f',Source_Num));
                     end          
                     Good(Source_Num) = 0; % mark it deleted
                     Goodndx(:,j) = []; % remove this column
                     % because of reverse sequencing, has no effect on the lower indices in the for loop
                  end
               end                  
            end
         end % for loop, testing adequacy of sources
      end % if this source is still good
      
   end % for all source combinations
   
   if(GUI.VERBOSE),
      close(hwait), end
   
   if GUI.EEG 
      % Detect EEG Reference 
      RefNdx = good_channel(Channel,[],'EEG REF');
      if isempty(RefNdx) % Average reference
         RefNdx = 0;
      end
   end
   
   if GUI.MEG
      ChanNdx = MEGndx;
   elseif GUI.EEG
      ChanNdx = EEGndx;    
   end
   
   HeadModel.Gain{iOrder} = NaN*ones(length(Channel),length(Goodndx(:)));
   
   if GUI.VERBOSE
      bst_message_window(sprintf('Decimating source grid to %.0f most significant sources',...
         length(find(Good)))); end
   
   HeadModel.Gain{iOrder}(ChanNdx,:) = G(:,Goodndx(:)); clear G Goodndx
   HeadModel.GridLoc{iOrder} = L(:,find(Good)); clear L Good
   
   % Now compute EEG gain matrix from MEG's source grid locations if MEG and EEG are available
   if GUI.MEG & GUI.EEG
      if GUI.VERBOSE
         bst_message_window(sprintf('Now creating EEG gain matrix for %.0f sensors and %.0f sources',...
            length(EEGndx),size(HeadModel.GridLoc{iOrder} ,2))); end
      
      
      if iOrder == 1 | iOrder == 4 % Current dipole source models
         G = feval(HeadModel.Function{EEGndx(1)},HeadModel.GridLoc{iOrder},Channel,HeadModel.Param,ORDER); % same for synchronous sources as well    
      else
         G = NaN*ones(size(HeadModel.Gain{iOrder}(EEGndx,:)));
      end
      
      HeadModel.Gain{iOrder}(EEGndx,:) = G;%G(EEGndx,:);     
      
   end
   
   clear G
   
   if 0, %GUI.VERBOSE
      bst_message_window({sprintf('Creating %s gain matrices for volume source scan -> DONE', SourceOrderString{iOrder}),' '})
   end
   
end % per source model order

% % for now, let's handle the synchronous sources somewhat crudely. I'm gonna
% %  just use the multipolar grid, and let rapmusic_gui sample pairs from that
% 
% for i = VALIDORDER
%     iOrder = i+2;
%     HeadModel.GridLoc{iOrder+3} = HeadModel.GridLoc{iOrder};
%     HeadModel.Gain{iOrder+3} = HeadModel.Gain{iOrder};
% end

return
