function h = plotgaussians(mu1,sigma1,mu2,sigma2,xrange,yrange,linespecs, ...
			   contourlines)
% h = plotgaussians(mu1,sigma1,mu2,sigma2,xrange,yrange,linespecs,contourlines)
% Plots contours of two 2D Gaussians
% mu1,sigma2 and mu2,sigma2 each specify a 2-dim Gaussian
% as well as the decision boundary between them, given by priors
% All other arguments are optional: 
% xrange=[fromx tox] specifies range of x values to plot (default -10 10)
% yrange=[fromy toy] specifies range of y values to plot (default same as x)
% linespecs is a cell array specifying line specifications for (in this order)
%   the decision boundary, the first Gaussian, the second Gaussian
% contourlines specifies the contour line values for the Gaussians

  if nargin<8
      linespecs = {'-k',':r','--b'};
  end
  if nargin<5
      xrange = [-5 5];
  end
  if nargin<6
      yrange = xrange;
  end
  if nargin<10
      contourlines = 5 ;
  end
  [x,y] = meshgrid(fillrange(xrange),fillrange(yrange));
  g1 = gaussian2d(x,y,mu1,sigma1);
  g2 = gaussian2d(x,y,mu2,sigma2);
  if length(contourlines)==1
      contourlines=0:(1.5/contourlines):1.5 ;
  end
  h = figure;
  hold;
  contour(x,y,g1*sqrt(det(sigma1)),contourlines,linespecs{2});
  contour(x,y,g2*sqrt(det(sigma2)),contourlines,linespecs{3});
  contour(x,y,g1-g2,[0 0],linespecs{1});
  axis equal;
function filledrange = fillrange(range)
  filledrange = range(1):((range(2)-range(1))/100):range(2);

function f = gaussian2d(x,y,mu,Sigma)
  % Calculated the 2-dimmensional Gaussian density at many points.
  % x and y must be matrices of the same shape. The result is also of the
  % same shape, and holds the corresponding density values.
  [n,m] = size(x) ;
  f = multigaussian([reshape(x,n*m,1) reshape(y,n*m,1)],reshape(mu,1,2),Sigma);
  f = reshape(f,n,m);
  
