function [gridlocs] = gridmaker(vertices,faces,grid_fname,Verbose)
%GRIDMAKER - 3D regular gridding of surface envelope
% function [gridlocs] = gridmaker(vertices,faces,grid_fname,Verbose)
% function [gridlocs] = bem_gridmaker(vertices,faces,bem_grid_mfname, OPTIONS)
% 3D regular gridding of the inside of surface envelope
% INPUTS
% vertices a cell containing the Nx3 set of tessellation vertices of the enveleppe
% faces    a cell containing the Mx3 array defining the triangles of the tessellated envelope from the set of vertices
% grid_fname character array where to store the grid locations  
%   Verbose  = 1 turns on Verbose mode
%  
%   OUTPUTS
%   gridlocs is a Rx3 array of grid point locations

%<autobegin> ---------------------- 08-Jun-2004 15:11:20 -----------------------
% --------- 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
%
% At Check-in: $Author: Mosher $  $Revision: 14 $  $Date: 6/08/04 1:38p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 28-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> ------------------------ 08-Jun-2004 15:11:20 -----------------------

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

SPACING = 5/1000; % in meters
NGridMax = 20000; % About twice the maximum number of grid points to achieve 
NGridMin = 10000; % About twice the minimum number of grid points 

SubjectFV.faces = faces{1};
SubjectFV.vertices = vertices{1};

if size(SubjectFV.faces,1)>500
    NFV = reducepatch(SubjectFV,500); % number of faces
else
    NFV = SubjectFV;
end

center = mean(NFV.vertices);
NFV.vertices = .9 * (NFV.vertices - repmat(center,size(NFV.vertices,1),1)); %.9 - don't put grid points too close to inner enveloppe
NFV.vertices = NFV.vertices + repmat(center,size(NFV.vertices,1),1);

GridBox = NFV.vertices'; % about half that for the vertices

mx = max(GridBox,[],2);
mn = min(GridBox,[],2);

if Verbose
    bst_message_window(...
        'Gridding the inner BEM envelope. . .'...
        )
end

L = zeros(3,NGridMax+1);

while size(L,2) > NGridMax | size(L,2) < NGridMin % Maximum number of grid points
    
    if size(L,2) > NGridMax
        SPACING = SPACING * 1.3;
    elseif size(L,2) < NGridMin
        SPACING = SPACING / 1.2;
    end
    
    vx = [mn(1):SPACING:mx(1)];
    vy = [mn(2):SPACING:mx(2)];
    vz = [mn(3):SPACING:mx(3)];

    [X,Y,Z] = meshgrid(vx,vy,vz);
    
    L = [X(:) Y(:) Z(:)]';
    
end

is = inpolyhd(L',NFV.vertices, NFV.faces); % find points within inner BEM envelope

Rq_bemgrid = L(:,find(is))'; % downselect

if Verbose
    
   bst_message_window(...
       'Gridding the inner BEM envelope -> DONE'...
       )
   
end


save(grid_fname,'Rq_bemgrid');

if nargout == 1
    gridlocs  = Rq_bemgrid;
end


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


