// Shrink the input array by removing unassigned elements
//
#macro ArrayShrink(A)
 #local N = dimension_size(A, 1);
 #local RN = 0;
 #for (I, 0, N-1)
  #ifdef (A[I]) #local RN = RN + 1; #end
 #end
 #if (RN = 0)
  #break // exit without changes if input array is empty
 #end
 #if (RN = N)
  #break // exit without changes if input array is fully assigned
 #end
 #local RI = 0;
 #local OutA = array[RN];
 #for (I, 0, N-1)
  #ifdef (A[I])
   #local OutA[RI] = A[I];
   #local RI = RI + 1;
  #end
 #end
 #declare A = OutA;
#end

// Create an array of 2D polygons of Voronoi cells based on the
//  input seed point array.
// Input point y coordinates are ignored.
// Output is an array the same size as the input array.
// Each element of the output array is a point array representing a cell polygon
//  on the x-z plane.
//
#macro VoronoiCells(Seeds, CellsOut)
 #local N = dimension_size(Seeds, 1);
 // Find extent of seed list
 #local MinX = Seeds[0].x; #local MinZ = Seeds[0].z;
 #local MaxX = Seeds[0].x; #local MaxZ = Seeds[0].z;
 #for (I, 1, N-1)
  #if (Seeds[I].x < MinX) #local MinX = Seeds[I].x; #end
  #if (Seeds[I].z < MinZ) #local MinZ = Seeds[I].z; #end
  #if (Seeds[I].x > MaxX) #local MaxX = Seeds[I].x; #end
  #if (Seeds[I].z > MaxZ) #local MaxZ = Seeds[I].z; #end
 #end
 // Make depth the size of extents to guarantee full cells
 #local D = max(MaxX-MinX, MaxZ-MinZ) / 2;
 // Make trace target unions for each seed point
 #local OtherCones = array[N];
 #for (I, 0, N-1)
  #local OtherCones[I] = union {
   #for (J, 0, N-1)
    #if (J != I)
     cone { <Seeds[J].x, -D, Seeds[J].z>, D, <Seeds[J].x, 0, Seeds[J].z>, 0 }
    #end
   #end
   plane { y, -D }
  }
 #end
 // Trace cells
 #local NC = 200;
 #local DA = 360 / NC;
 #local Cells = array[N];
 #for (I, 0, N-1)
  #local Cells[I] = array[NC];
  #for (J, 0, NC-1)
   #local V = vrotate(<1,-1,0>, y*DA*J);
   #local P = trace(OtherCones[I], <Seeds[I].x, 0, Seeds[I].z>, V);
   #local P = <P.x, 0, P.z>;
   #local Cells[I][J] = P;
  #end
 #end
 // Remove collinear points
 #for (I, 0, N-1)
  #local NewCell = array[NC];
  #local NewI = 0;
  #for (J, 0, NC-1)
   #local JB = (J = 0 ? NC-1 : J-1);
   #local JA = (J = NC-1 ? 0 : J+1);
   #local VB = vnormalize(Cells[I][J] - Cells[I][JB]);
   #local VA = vnormalize(Cells[I][JA] - Cells[I][J]);
   #if (abs(vdot(VB, VA)-1) > 1e-4)
    #local NewCell[NewI] = Cells[I][J];
    #local NewI = NewI + 1;
   #end
  #end
  ArrayShrink(NewCell)
  #local Cells[I] = NewCell;
 #end
 // Return array of cell polygons
 #declare CellsOut = Cells;
#end