function [REL_AREAS] = voronoiGenerator2d(M,p);
% MYVORONOI(M,p) ----
% This function will take standardized point configurations, and make voronoi 
% tessellations of them.  The edges are bounded artificially, by putting an octagon of
% points around the field, far outside of it.  These will absord the infinited edge
% cells, leaving the points of interest with finite areas.  This will not 
% come out perfectly or consistently shaped.  So, RELATIVE "domains of danger" will 
% be found, since the tessellation will have different total areas for each case.  At 
% least the areas are finite.  
%
% p is an optional input; if it is 'y', then the voronoi diagram will be plotted
%  after calculation.  Pretty to look at, but the plot command is a pain.
if nargin == 1
    plot_req = 0; 
elseif strcmp(p,'y') ==1
    plot_req = 1;
end

% First, some parameters of the field.  The 'width' is the actual half-length of the 
% field dimensions; and L will be the distance out that the artificial
%  edges occur. 
% W = sqrt(9*100*pi);
% L = W*4;
L = 1;
OB = (1/sqrt(2)) * L;

% The boundary is made by placing an octogon around the field, far outside of it.
% BOUND is a matrix containing the coordinates of the octogon vertices.
p = size(M)

if p(2) == 4
   BOUND = [0 0 -L 0; 0 0 L 0; 0 -L 0 0; 0 L 0 0; 0 -OB -OB 0; 0 OB OB 0; 0 OB -OB 0;...
         0 -OB OB 0];
else
   BOUND = [0 0 -L ; 0 0 L ; 0 -L 0 ; 0 L 0 ; 0 -OB -OB ; 0 OB OB ; 0 OB -OB ;...
         0 -OB OB];
end

% find the center of mass, using the centerOfMass2d function
    
center = [];
center = centerOfMass2d(M);
centerX = center(1,1);
centerY = center(1,2);

% find the greatest distance between any two points, which is going to
% be our group diameter
    
currentMax = 0;
j=0;

for i=1:(size(M,1))
	for j=1:(size(M,1))
		if dist2d(M(i,2), M(i,3),  M(j,2), M(j,3))>currentMax
                currentMax = dist2d(M(i,2), M(i,3), M(j,2), M(j,3));
    	end
  	end
end
    
diameter = currentMax;
    
% Now we blow the octagon up to make sure that it's larger than our diameter

BOUND = BOUND*diameter*1.5; 

% Now we recenter it at the center of mass.

for i = 1:size(BOUND,1)
 	BOUND(i,2)=BOUND(i,2)+centerX;
   BOUND(i,3)=BOUND(i,3)+centerY;
end

    
    

% The edge points are added to the input array M, then the x and y coordinates 
% extracted from columns 2 and 3, since this 2-column format is needed
% for the 'voronoin' function.
N0 = [M; BOUND];
N = [N0(:,2) N0(:,3)];


% Now, the voronoi calculations can be made.  The output V will be a matrix containing
% the coordinates of all vertices for the polygons; C is a cell, each element 
% (rowwise) being a matrix containing the indices into V of the necessary vertices 
% for the polygon around that point. 
[V C]  = voronoin(N);
    

% Then, for each point i, the polygon is constructed from the voronoi data in C and
% V, and the area computed by convhulln.  K will just be a repeat of the vertices, 
% but A will be the area and therefore will be saved.  By basing the no. of reps on 
% size(M), it will ignore the tessellations around the edge points (i.e., the 
% infinite throwaways).

AREAS = [];                             % Initialize matrix for area storage.
for i = 1:size(M,1)                     % Loop over each point.
    
    X = V(C{i},:);                      % All the necessary vertices are put into 
                                        %  matrix X.
    X
    [K A] = convhull(X(:,1),X(:,2));               % Area is computed for the polygon.
    
    AREAS = [AREAS; A];                 % Save the areas into the matrix.
end

% Areas are added up, and the 'RELative AREAS' matrix initialized, for storage of, 
% you guessed it, the relative areas of each polygon, for each point.
TOTAL_A = sum(AREAS);


REL_AREAS = [];

% Calculate relative areas, and store them.
for j = 1:size(AREAS,1)                 % Loop through each polygon.
    
    relarea = AREAS(j,1) / TOTAL_A;     % Calculate rel a.
    
    REL_AREAS = [REL_AREAS; relarea];   % Store each in the matrix.
end

% Values are multiplied by 100 (i.e. put in % form) to make them easier 
% to read.
REL_AREAS = REL_AREAS * 100;


hold on
% If a plot was requested, the following commands will be enabled.
if plot_req == 1
  p = voronoi(N0(:,2),N0(:,3));     % This command actually generates a 
                                    %  voronoi diagram, named 'p'.
                                    %axis([-60,60,-60,60])         % Axes are formatted.
                                    %!!!!NOTE Right now I'm letting Matlab 
                                    %!!!!     pick its own boundaries, but
                                    %!!!!     we should probably help it a little bit
                                    %!!!! 	 later on.  UPDAATED 22JUNE2005 Jim Work
                                    % Set the color and sizes for the actual lines and dots on the graph.
                                    
  set(p,'markersize',4,'markeredgecolor','k','markerfacecolor','k', 'linewidth', 2)
 
  % Set the axis and background colors and properties. 
  set(gca,'color','w','xcolor','k','ycolor','k','tickdir','out','fontweight','bold',...
      'fontsize',12);
  set(gcf,'color','none');

  
  for h = p
    set(h,'color','k','markersize',8);
  end
    
  plot(M(:,2),M(:,3),'o','color','k','markerfacecolor','k','markersize',4)
  plot(BOUND(:,2), BOUND(:,3), 'r*') % NOTE: To mask the bound markers, change r to w on this line.
end


hold off