function GaussianRegionsDisplayUtility(meanstruct,Cstruct,classpriors,data)
%function GaussianRegionsDisplayUtility(meanstruct,Cstruct)
% Displays 2-D acceptance regions for the classes specified.
% meanstruct:  [mu1 mu2 ... muN] is a matrix of mean vectors
% Cstruct: 2 x 2 x N array  of covariance matrices --
%  C1 = Cstruct(:,:,1) ...   CN = C(:,:,N)
% classpriors: 1 x N vector of classprobabilities
% data:  structure array of datamatrices --  data(1).pts .... data(N).pts
%  data(1).pts is 2 x Mc, where Mc is the number of points for the cth class
%
% Utility function to help students visualize the class boundaries for 2-D Gaussian
% decision problems
%  For CSCI 5521, Fall 2003, University of Minnesota
%  Last modification Fall 2007
% author-  P.R.  Schrater

if length(classpriors)~=size(meanstruct,2) | length(classpriors)~=size(Cstruct,3) | size(meanstruct,2)~=size(Cstruct,3),
	error('Argument dimensions do not match');
end
maxvars = [0 0];
for j=1:size(Cstruct,3)  % find  bounds for sampling
 	[U,S,V]=svd(Cstruct(:,:,j));
	maxvars = max(maxvars,max(S));
end
upperbounds = max(meanstruct')+2*maxvars;
lowerbounds = min(meanstruct')-2*maxvars;

x = lowerbounds(1):(upperbounds(1)-lowerbounds(1))/255:upperbounds(1);
y = lowerbounds(2):(upperbounds(2)-lowerbounds(2))/255:upperbounds(2);

[X,Y]=meshgrid(x,y);

for j=1:length(classpriors),
	z1 = [X(:)-meanstruct(1,j) Y(:)-meanstruct(2,j) ]';
	pz(j,:) = classpriors(j)/(2*pi*sqrt(det(Cstruct(:,:,j))))*exp(-0.5*sum(z1.*(inv(Cstruct(:,:,j))*z1)));
end;
 [r,index]=sort(pz);  maxclass = index(end,:); 
 % convert class indices into colors
 im = reshape(220*(maxclass)/length(classpriors),size(X));
 h = hsv(256);
 colormap(h)
image(x,y,im); hold on;
if nargin>3,
	for j=1:length(classpriors),
		scatter(data(j).pts(1,:),data(j).pts(2,:),4*ones(size(data(j).pts(2,:))),h(j*round(220/length(classpriors))+35,:));
	end
end

hold off