function [vol, sens] = headmodelplot(cfg, data);

% HEADMODELPLOT makes a 3D visualisation of the volume conductor model
% and optionally of the gradiometer positions and headshape. It can
% be used for example to check CTF multiple-sphere head models.
%
% Use as
%   headmodelplot(cfg)
%   headmodelplot(cfg, data)
%
% The configuration should specify the head model
%   cfg.vol         = volume conductor definition
%   cfg.hdmfile     = filename
%
% For MEG, the gradiometer positions can be stored in the data.
% Otherwise, the gradiometers should be specified using either
%   cfg.gradfile    = filename
%   cfg.grad        = gradiometer defintion
%
% For EEG, the electrode positions can be stored in the data.
% Otherwise, the electrodes should be specified using either
%   cfg.elecfile    = filename
%   cfg.elec        = electrode defintion
%
% Other options are
%   cfg.channel          = cell-array, see CHANNELSELECTION
%   cfg.plotheadsurface  = 'yes' or 'no', is constructed from head model
%   cfg.headshape        = name of CTF *.shape file
%   cfg.plotspheres      = 'yes' or 'no'
%   cfg.plotspherecenter = 'yes' or 'no'
%   cfg.plotsensors      = 'yes' or 'no' plot electrodes or gradiometers
%   cfg.plotcoil         = 'yes' or 'no' plot all gradiometer coils
%   cfg.plotlines        = 'yes' or 'no' plot lines from sensor to head surface
%   cfg.spheremesh       = number of vertices for spheres, either 42, 162 or 642

% Copyright (C) 2004, Robert Oostenveld
%
% $Log: headmodelplot.m,v $
% Revision 1.7  2005/12/14 10:44:23  roboos
% implemented support for EEG, spherical and BEM
% implemented support for single sphere in  case of MEG
% fixed some incompatibilities in grad and vol structure for MEG
% removed cfg.plotheadshape, just look whether isempty(cfg.headshape)
%
% Revision 1.6  2005/11/24 16:26:07  roboos
% updated multisphere selection to reflect the current change in multisphere vol structure (no labels any more)
%
% Revision 1.5  2005/05/17 17:50:37  roboos
% changed all "if" occurences of & and | into && and ||
% this makes the code more compatible with Octave and also seems to be in closer correspondence with Matlab documentation on shortcircuited evaluation of sequential boolean constructs
%
% Revision 1.4  2004/08/19 07:22:49  roboos
% added cfg.plotspherecenter option
%
% Revision 1.3  2004/08/19 07:19:55  roboos
% added cfg.plotcoil option
%
% Revision 1.2  2004/08/06 12:37:25  roboos
% fixed minor bugs
%
% Revision 1.1  2004/08/06 08:55:23  roboos
% new implementation
%

% these are suitable RGB colors
skin   = [255 213 119]/255;
brain  = [202 100 100]/255;
cortex = [255 213 119]/255;
surface_facecolor = skin;
surface_edgecolor = 'none';
surface_facealpha = 0.7;

if ~isfield(cfg, 'vol') && ~isfield(cfg, 'hdmfile')
  cfg.vol = [];
end

if nargin<2
  data = [];
end

% set the defaults that apply both to EEG and MEG
if ~isfield(cfg, 'spheremesh'),       cfg.spheremesh = 642;          end
if ~isfield(cfg, 'plotsensors'),      cfg.plotsensors = 'yes';       end
if ~isfield(cfg, 'plotheadsurface'),  cfg.plotheadsurface = 'yes';   end

% extract/read the gradiometer and volume conductor
[vol, sens, cfg] = prepare_vol_sens(cfg, data);

% determine the type of input data
ismeg          = isfield(sens, 'ori');
iseeg          = (~ismeg);
isbem          = isfield(vol, 'bnd');
issphere       = isfield(vol, 'r');
ismultisphere  = isfield(vol, 'r') && length(vol.r)>4;
issinglesphere = isfield(vol, 'r') && length(vol.r)==1;

if ismeg
  % sensors describe MEG data, set the corresponding defaults
  if ~isfield(cfg, 'plotspheres'),      cfg.plotspheres = 'no';        end
  if ~isfield(cfg, 'plotspherecenter'), cfg.plotspherecenter = 'no';   end
  if ~isfield(cfg, 'plotlines'),        cfg.plotlines = 'yes';         end
  if ~isfield(cfg, 'plotcoil'),         cfg.plotcoil = 'no';           end
  if ~isfield(cfg, 'headshape'),        cfg.headshape = [];            end
else
  % sensors describe EEG data, set the corresponding defaults
  if ~isfield(cfg, 'plotspheres'),      cfg.plotspheres = 'no';        end
  if ~isfield(cfg, 'plotspherecenter'), cfg.plotspherecenter = 'no';   end
  if ~isfield(cfg, 'plotlines'),        cfg.plotlines = 'yes';         end
end

if ismeg && strcmp(megsystem(sens), 'ctf151')
  Nsensors = 151;
else
  Nsensors = length(sens.pnt);
end

if issphere
  % determine the number of spheres in the volume model
  Nspheres = length(vol.r);
  fprintf('spherical head model with %d spheres\n', Nspheres);
end

clf
hold on
axis equal
axis vis3d
if ismeg
  axis off
else
  axis on
  grid on
  xlabel('x')
  ylabel('y')
  zlabel('z')
end

if iseeg
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % plotting for EEG
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if strcmp(cfg.plotsensors, 'yes')
    plot3(sens.pnt(1:Nsensors,1), sens.pnt(1:Nsensors,2), sens.pnt(1:Nsensors,3), 'b*');
  end % plotsensors

  if strcmp(cfg.plotheadsurface, 'yes')  && ~isempty(vol)
    [pnt, tri] = headsurface(vol, sens);
    h = patch('Faces', tri, 'Vertices', pnt);
    set(h, 'facecolor', surface_facecolor);
    set(h, 'edgecolor', surface_edgecolor);
    set(h, 'facealpha', surface_facealpha);
  end % plotheadsurface

  if strcmp(cfg.plotlines, 'yes') && ~isempty(vol)
    if isbem
      % project the electrodes on the skin surface, on the nearest triangle
      [el] = project_elec(sens.pnt, vol.bnd(vol.skin).pnt, vol.bnd(vol.skin).tri);
      % this returns [tri, la, mu]
      tri = el(:,1);
      la  = el(:,2);
      mu  = el(:,3);
      for i=1:Nsensors
        v1  = vol.bnd(vol.skin).pnt(vol.bnd(vol.skin).tri(tri(i),1),:);
        v2  = vol.bnd(vol.skin).pnt(vol.bnd(vol.skin).tri(tri(i),2),:);
        v3  = vol.bnd(vol.skin).pnt(vol.bnd(vol.skin).tri(tri(i),3),:);
        prj(i,:) = routlm(v1, v2, v3, la(i), mu(i));
      end
    elseif issphere
      % project the electrodes onto the sphere surface, towards the origin
      nrm = sqrt(sum(sens.pnt.^2,2));
      prj = vol.r(vol.skin) * sens.pnt ./ [nrm nrm nrm];
    else
      % in case that no known volume conductor is specified
      prj = sens.pnt;
    end
    x = [sens.pnt(:,1) prj(:,1)]';
    y = [sens.pnt(:,2) prj(:,2)]';
    z = [sens.pnt(:,3) prj(:,3)]';
    line(x, y, z, 'color', 'm');
  end % plotlines

elseif ismeg
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % plotting for MEG
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if strcmp(cfg.plotsensors, 'yes')
    % plot only the bottom coil
    plot3(sens.pnt(1:Nsensors,1), sens.pnt(1:Nsensors,2), sens.pnt(1:Nsensors,3), 'b*');
  end % plotsensors

  if strcmp(cfg.plotcoil, 'yes')
    pnt = sens.pnt;
    ori = sens.ori;
    x = [pnt(:,1) pnt(:,1)+ori(:,1)]';
    y = [pnt(:,2) pnt(:,2)+ori(:,2)]';
    z = [pnt(:,3) pnt(:,3)+ori(:,3)]';
    line(x, y, z, 'color', 'g');
    plot3(pnt(:,1), pnt(:,2), pnt(:,3), 'g.');
  end % plotsensors

  if ~isempty(cfg.headshape)
    [shape] = read_ctf_shape(cfg.headshape);
    plot3(shape.pnt(:,1), shape.pnt(:,2), shape.pnt(:,3), 'r.');
  end % plotheadshape

  if strcmp(cfg.plotspheres, 'yes') && ~isempty(vol)
    % create a triangulated unit sphere
    if cfg.spheremesh==42
      [pnt0, tri] = icosahedron42;
    elseif cfg.spheremesh==162
      [pnt0, tri] = icosahedron162;
    elseif cfg.spheremesh==642
      [pnt0, tri] = icosahedron642;
    end
    Nvertices = size(pnt0,1);
    for i=1:Nspheres
      % scale and shift the unit sphere to the proper location
      pnt = pnt0*vol.r(i) + repmat(vol.o(i,:),Nvertices,1);
      if Nspheres>4
        distance = sqrt(sum((pnt - repmat(sens.pnt(i,:),Nvertices,1)).^2, 2));
        distance = (distance-min(distance))/range(distance);
      else
        distance = zeros(Nvertices, 1);
      end
      h = patch('Faces', tri, 'Vertices', pnt);
      % set(h, 'FaceVertexCData', 0.5*ones(length(distance),30));
      % set(h, 'FaceVertexCData', [0 0 1]);
      set(h, 'FaceVertexAlphaData', (1-distance).*(distance<0.1));
      % set(h, 'AlphaDataMapping', 'direct');
      set(h, 'edgealpha', 'interp')
      % set(h, 'edgealpha', 1)
      set(h, 'facealpha', 'interp')
      % set(h, 'facealpha', 0)
      set(h, 'facecolor', 'none');
      % set(h, 'linestyle', 'none');
    end
  end % plotspheres

  if strcmp(cfg.plotspherecenter, 'yes') && ~isempty(vol)
    plot3(vol.o(:,1), vol.o(:,2), vol.o(:,3), 'k.');
  end % plotspherecenter

  if strcmp(cfg.plotheadsurface, 'yes') && ~isempty(vol)
    % estimate the head surface from the spheres and gradiometers
    [pnt, tri] = headsurface(vol, sens);
    h = patch('Faces', tri, 'Vertices', pnt);
    set(h, 'facecolor', surface_facecolor);
    set(h, 'edgecolor', surface_edgecolor);
    set(h, 'facealpha', surface_facealpha);
  end % plotheadsurface

  if strcmp(cfg.plotlines, 'yes') && ismultisphere
    % first determine the indices of the relevant gradiometers.
    % [sel_g, sel_v] = match_str(sens.label, vol.label);
    sel_g = 1:Nsensors;
    sel_v = 1:Nsensors;
    % create head-surface points from multisphere head-model.
    dir     = sens.pnt(sel_g,:) - vol.o(sel_v,:);
    dist    = sqrt(sum(dir.*dir,2));
    pnt0    = repmat((vol.r(sel_v)./dist),1,3).*dir + vol.o(sel_v,:);
    pnt1    = sens.pnt(sel_g,:);
    f = gcf;
    figure
    hist(sqrt(sum((pnt0-pnt1).^2, 2)));
    xlabel('distance (a.u.)');
    title('histogram of distances from the bottom coil to the head surface');
    figure(f);
    x = [pnt0(:,1) pnt1(:,1)]';
    y = [pnt0(:,2) pnt1(:,2)]';
    z = [pnt0(:,3) pnt1(:,3)]';
    line(x, y, z, 'color', 'm');
  end % plotlines

end % iseeg or ismeg

lighting gouraud
camlight
hold off
