function [REL_VOLS] = voronoiGenerator3d(M);
% myvoronoi3d(M) ----  Calculates Voronoi volumes for an array of points
% in 3D.  Analogous to the 'myvoronoi' function, with an extra column in 
% the input for z-coordinates.  Here also, the points-of-interest are 
% bounded far outside the field-of-interest to avoid having infinite volumes.
% The vertices of a dodecahedron form the boundary.
%
% The input matrix 'M' should have individual numbers in the 1st column, then
% x y and z coordinates in the 2nd, 3rd, and 4th columns respectively.
%
% Output are the relative volumes (in a vector) of the polyhedra for each 
% point, i.e. the percent-area each occupies.


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

% Next, the boundary points will be placed, far outside of the space of interest.  
% The vertices of a dodecahedron are used, 20 points total, and are placed in the 
% 'BOUND' matrix.
BOUND = [148.3282 -388.3282   0;...
 -148.3282 -388.3282         0;...
 -240.0000 -240.0000  240.0000;...
         0 -148.3282  388.3282;...
  240.0000 -240.0000  240.0000;...
         0  148.3282 -388.3282;...
 -240.0000  240.0000 -240.0000;...
 -148.3282  388.3282         0;...
  148.3282  388.3282         0;...
  240.0000  240.0000 -240.0000;...
         0 -148.3282 -388.3282;...
 -240.0000 -240.0000 -240.0000;...
 -388.3282         0 -148.3282;...
 -388.3282         0  148.3282;...
 -240.0000  240.0000  240.0000;...
         0  148.3282  388.3282;...
  240.0000  240.0000  240.0000;...
  388.3282         0  148.3282;...
  388.3282         0 -148.3282;...
  240.0000 -240.0000 -240.0000];

% A fourth column is added to 'BOUND', so that the input matrix M may be
% attached correctly.  
BOUND = [zeros(20,1) BOUND];



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


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

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

VOL = [];                             % 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.
    
    [K A] = convhulln(X);               % Volume is computed for the polyhedron.
    
    VOL = [VOL; A];                 % Save the volumes into the matrix.
end

% Volumes are added up, as TOTAL_V,
% and the 'RELative VOLumes' matrix initialized, for storage of
% the relative areas of each polyhedron, for each point.
TOTAL_V = sum(VOL);

REL_VOLS = [];

% Calculate relative volumes, and store them.
for j = 1:size(VOL,1)                 % Loop through each polyhedron.
    
    relarea = VOL(j,1) / TOTAL_V;     % Calculate rel a.
    
    REL_VOLS = [REL_VOLS; relarea];   % Store each in the matrix.
end


REL_VOLS = REL_VOLS * 100;