function Sphere = overlapping_sphere(Channel,SubjectFV, VERBOSE, GRAPHICS)
%OVERLAPPING_SPHERE - Create overlapping spheres for E/MEG channels.
% function Sphere = overlapping_sphere(Channel,SubjectFV, VERBOSE, GRAPHICS)
% 
% Channel is an array of structures.
% Channel.Loc, a matrix of locations of the sensor coil centers, one column per coil
% Channel.Orient, corresponding matrix of orientations, null in the case of EEG
% Channel.Weight, relative weights to each coil. Used here in the least-squares (default is 1).
% Channel.Type, 'MEG','EEG' are the only supported types
% Channel.Comment, unused here.
%
% SubjectFV.vertices contains the vertices information (3 x Nverts)
% SubjectFV.faces contains the corresponding faces information (Nfaces x 3)
%  The faces information is used for determining the surface normals
%  and for visualization.
%
% VERBOSE = 1 (optional), text info.
% GRAPHICS = 1 (optional), graphics displayed
%
% Sphere is an array of structures
% For each Channel, returns 
%  Sphere.Center, a 3 x 1 center location and 
%  Sphere.Radius, the corresponding radius.
%  Sphere.Weight is a weighting vector for coloring the vertices to indicate
%  the relative surface weights.
%  Sphere.Approx is the identified surface weights on the sphere.
%
%  Sphere.Center and Sphere.Radius should be transferred to fields in the 
%   structure HeadModel.Param

%<autobegin> ---------------------- 16-Jun-2004 16:50:21 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Forward Modeling
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bst_message_window.m
%   toolbox\colnorm.m
%   toolbox\overlapping_sphere_fmins.m
%   toolbox\plot3d.m
%   toolbox\source_grid.m
%   toolbox\tesselate.m
%   toolbox\view_surface.m
%   toolbox\weighting_scalar.m
%
% At Check-in: $Author: Mosher $  $Revision: 19 $  $Date: 6/16/04 3:16p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 16-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> ------------------------ 16-Jun-2004 16:50:21 -----------------------


% /---Script Author--------------------------------------\
% |                                                      |
% | *** John C. Mosher, Ph.D.                            |
% |  Biophysics Group                                    |
% \------------------------------------------------------/

% History ---------------------------------------------------- 
% Date of creation: Janurary 2001
% Date of modification: March 2002                  
%
% March 22, 2002 : SB Altered: Script header / response to VERBOSE mode
%                  / replaced FMINS by FMINSEARCH
% JCM 08-Sep-2003  Commenting
% JCM 18-Nov-2003 Replacing .Vertices and .Faces with .vertices and .faces 
%                 in the structure SubjectFV to be
%                 consistent with the bst_headmodeler call.
% SB  15-Jun-2004 Minor graphics enhancements
% ------------------------------------------------------------


% Calculate the Normals

if ~exist('VERBOSE','var') | isempty(VERBOSE),
   VERBOSE = 0; % default is silent running
end

if ~exist('GRAPHICS','var') | isempty(GRAPHICS),
   GRAPHICS = 0; % default is no display
end

%figure(windfind('Junk'))
% Matlab has the unusual spec of vertices being n x 3, not 3 x n, in patch.
hp = patch('faces',SubjectFV.faces,'vertices',SubjectFV.vertices','Visible','off');
Normals = get(hp,'VertexNormals');
[ignore,Normals] = colnorm(Normals'); % transpose and unit norm of each vertex
delete(hp)
%close(gcf)

Nc = length(Channel); % how many channels to process
% allocate
[Sphere(1:Nc)] = deal(struct('Center',[],'Radius',[],'Weight',[],'Approx',[]));

if VERBOSE
    hwait = waitbar(0,sprintf('Fitting Spheres to %.0f Channels',Nc));
    set(hwait,'units','norm'); % set units to normalized
    hwait_pos = get(hwait,'position');
    set(hwait,'position',[hwait_pos(1) 0.1 hwait_pos(3:4)]); % shift to bottom away from other images
    drawnow
end

for i = 1:Nc, % for each channel
    % Calculate the surface weight for the true shape
    if isempty(Channel(i).Weight)
        Channel(i).Weight = 1;    
    end
    
    Sphere(i).Weight = weighting_scalar(Channel(i).Loc,Channel(i).Orient,...
        Channel(i).Weight,SubjectFV.vertices, Normals);
    
    % initialize the guess of a sphere center as the weighted average of the
    % vertices
    %  tempCenter = mean(SubjectFV.vertices.*[[1;1;1]*Sphere(i).Weight],2); % weighted avg
    tempCenter = mean(SubjectFV.vertices,2);
    
    % distance from sensor to center is
    tempD = colnorm(Channel(i).Loc(:,1) - tempCenter);
    
    if(0)
        % initialize radius as average distance to center
        tempR = mean(colnorm(...
            SubjectFV.vertices - tempCenter(:,ones(1,size(SubjectFV.vertices,2)))));
    else
        % 11/5/99 above was nice idea, but sometimes yielded radii larger than sensor distance,
        %  due to neck region. So keep it easy
        tempR = min(.9*tempD,abs(tempD - 0.03)); % come in 3 cm, handle pathological cases
    end
    
    % grid something comparable
    testL = source_grid(tempD,tempR,1/1000,.3)'; % sparse sampling
    testL = [testL [testL(1:2,:);-testL(3,:)]]; % upper and lower hemispheres
    
    testL = testL + tempCenter(:,ones(1,size(testL,2)));
    
    % scan over this grid for weighting scalars, find the best match
    
    BestGridPt = 1; % assume first
    best_err = overlapping_sphere_fmins([testL(:,1);tempR],...
        Sphere(i).Weight,Channel(i).Loc,Channel(i).Orient,Channel(i).Weight,SubjectFV.vertices);
    for iGrid = 2:size(testL,2), % for each additional grid point
        test_err = overlapping_sphere_fmins([testL(:,iGrid);tempR],...
            Sphere(i).Weight,Channel(i).Loc,Channel(i).Orient,Channel(i).Weight,SubjectFV.vertices);
        if(test_err < best_err),
            BestGridPt = iGrid;
        end
    end
    
    % set start of search here.
    tempCenter = testL(:,iGrid);
    tempR = mean(colnorm(...
        SubjectFV.vertices - tempCenter(:,ones(1,size(SubjectFV.vertices,2)))));

    if VERBOSE 
        bst_message_window({...
                sprintf('Channel # %d / %d',i,length(Channel)),...    
                sprintf('Initial location is %2.2f %2.2f %2.2f (cm)',100*tempCenter),...
                sprintf('Initial radius is %2.2f (cm)',100*tempR)})
    end
        
    
    X = fminsearch('overlapping_sphere_fmins',[tempCenter;tempR],[],...
        Sphere(i).Weight,Channel(i).Loc,Channel(i).Orient,Channel(i).Weight,SubjectFV.vertices);
    Sphere(i).Center = X(1:3); % the location
    
    [err,sc,Radius] = overlapping_sphere_fmins(X,...
        Sphere(i).Weight,Channel(i).Loc,Channel(i).Orient,Channel(i).Weight,SubjectFV.vertices);
    Sphere(i).Radius = Radius; % the radius, mapped back depending on the Case Switch. 
    Sphere(i).Approx = sc;
    
    
    if VERBOSE,
       
       if GRAPHICS & i == 1,
          fig_os = figure;
          set(fig_os,'Name','Overlapping-Sphere Fitting')
          r = ceil(sqrt(length(Channel))); % Number of rows for following subplot
          c = ceil(sqrt(length(Channel)));
          FigName = get(fig_os,'Name');
       end
       
       bst_message_window({...
             sprintf('Final location is %2.2f %2.2f %2.2f (cm)',100*Sphere(i).Center),...
             sprintf('Final radius is %2.2f (cm)',100*Sphere(i).Radius),' '});
       
        
       if GRAPHICS,
          figure(fig_os)
          axx = subplot(r,c,i);
          if i == 1
              [hf,hs,hl] = view_surface(FigName,...
                  SubjectFV.faces,SubjectFV.vertices,abs(Sphere(i).Weight));
              set(hs,'facealpha',.8)
              set(hs,'edgecolor','flat');
              camlight, lighting gouraud    
          else % Spare time by just copying the basic scalp objects
              C = copyobj([hs,hl],axx);
              axes(axx)
              axis off, 
              hs = C(1);
          end
          hold on
          [ignore,ignore,ignore,SVerts,SFaces] = tesselate(Sphere(i).Radius,.15*Sphere(i).Radius,'full');
          SVerts = SVerts'; % columns are vertices
          [ignore,SNormals] = colnorm(SVerts);
          SVerts = SVerts + Sphere(i).Center*ones(1,size(SVerts,2));  
          % Channel graph object
          p=patch('vertices',SVerts','faces',SFaces','facecolor',[1 .8 .8],'edgecolor',[1 .5 .5]);
          plot3d(Channel(i).Loc','g+')
          set(p,'Facealpha',.6,'edgecolor','none')
          
          axis vis3d tight equal
          view(-90,0);
          drawnow
          
          figure(hwait) % bring to foreground
       end
       
       waitbar(i/Nc,hwait); % update waitbar
       drawnow
       
    end
    
end

if(VERBOSE),
    close(hwait);
end

return

