function [h, h3D] = plotmixture2D(X, evalfunc, compParams, varargin)
  
  showclassifications = 1;
  show3D = 1;
  PointSize = 1;
  PointStyle = '.';
  CompLinespec = 'b-.';
  CompContourLines = 5;
  DispPrior = 1;
  DensityLinespec = 'r-';
  setcallervars(varargin);
  if length(CompContourLines)==1
    CompContourLines=0:(1.5/CompContourLines):1.5;
  end
  
  mm = [min(X) ; max(X)] ;
  mm = mm + [-0.05 ; 0.05] * (mm(2,:)-mm(1,:)) ;
  x1range = mm(:,1);
  x2range = mm(:,2);
  [x1,x2] = meshgrid(fillrange(x1range),fillrange(x2range));
  meshshape = size(x1);
  x1x2 = [ reshape(x1,[],1) reshape(x2,[],1) ] ;

  % Evaluate the actuall mixture density
  p = feval(evalfunc,x1x2,compParams);
  p = reshape(p,meshshape);
  % And the most likely component of each point of X
  if showclassifications
    [dummy, c] = feval(evalfunc, X, compParams);
  else
    c=1;
  end
  
  if show3D
    h3D = figure;
    mesh(x1,x2,p);
  end
  h = figure;
  axes;
  axis([x1range ; x2range]);
  hold;
  % Plot the points of X
  scatter(X(:,1),X(:,2),PointSize,c,PointStyle);
  % Plot each component
  for i=1:length(compParams)
    g = multigaussian(x1x2,compParams(i).mean,compParams(i).cov) ;
    g = g*sqrt(det(compParams(i).cov)); % To get CDF instead of density
    g = reshape(g,meshshape);
    if CompLinespec
      contour(x1,x2,g,CompContourLines,CompLinespec);
    end
    if DispPrior
      text(compParams(i).mean(1),compParams(i).mean(2),...
	   num2str(compParams(i).prior),...
	   'HorizontalAlignment','center');
    end
  end
  % Plot the density contours
  if DensityLinespec
    contour(x1,x2,p,DensityLinespec);
  end
  

function filledrange = fillrange(range)
  filledrange = range(1):((range(2)-range(1))/100):range(2);
      