; poly_voronoi.pro   written Fall 2011 by Ken Desmond
;
; some comments added Aug 21-25, 2025 by Eric Weeks

; FUNCTION poly_voronoi, x, w, ri, ob, tri, threshold = threshold, Box = Box, Boundary = Boundary, area = area
;
; Box = [minx, maxx, miny, maxy]
;       default:  taken from the datao
; Boundary = two letters, R for rigid or P for periodic
;       default:  ['R', 'R']
; x = the positions in (x,y)
; w = the weights, typically the particle radii
; threshold = any Delaunay triangle that has an edge
;     exceeding this threshold will be deleted.  Presumably
;     these are triangles that are heading to infinity.  The
;     default is 8 times the largest radius, or 8 times
;     the typical inter-particle spacing, whichever is greater.
;
; RETURNS:
;  function itself returns the vertices of the Voronoi polygons
;
;  ri = reverse indices, for drawing the polygons
;  ob = 1 if a particle is on the edge, otherwise 0
;  tri = vertices of Delaunay triangles
;  area = area:  returns the areas of the Voronoi polygons
;   
; EXAMPLE:
;
; Let's consider an array 'a' that is 3xN with columns (x,y,r).
; We want to calculate the radical Voronoi polygons and plot
; them.
; 
; pv = poly_voronoi(a(0:1,*),a(2,*),ri,ob,tri)
;
; plot,a(0,*),a(1,*),/nodata,/iso
; na = n_elements(a(0,*))
; cx=cos(findgen(101)/100.*2.0*3.14159)
; sx=sin(findgen(101)/100.*2.0*3.14159)
; ;draw the circles:
; for i=0,na-1 do oplot,a(0,i)+cx*r[i],a(1,i)+sx*r[i]
; ;draw the polygons:
; for i=0,na-1 do begin & $
;   idx = ri[ri[i]:ri[i+1]-1] & $
;   idx = [idx,idx[0]] & $
;   oplot,pv(0,idx),pv(1,idx) & $
; endfor


; put a helper function first -- Ken wrote this as
; stand-alone, Eric is incorporating it here

FUNCTION kd_length, a, dim

s = size(a)

IF s[0] EQ 0 THEN return, 0

IF n_elements(dim) EQ 0 THEN BEGIN

   return, s[1]

ENDIF ELSE BEGIN

   IF dim GT s[0] THEN BEGIN
      return, 0
   ENDIF ELSE BEGIN
      return, s[dim]
   ENDELSE

ENDELSE

END

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

FUNCTION poly_voronoi, x, w, ri, ob, tri, threshold = threshold, Box = Box, Boundary = Boundary, area = area

w = reform(w)
n = kd_length(x, 2)
ob = make_array(n, value=0)

IF ~keyword_set(Box) THEN BEGIN 
   Box = [min(x[0, *]), max(x[0, *]), min(x[1, *]), max(x[1, *])]
   Boxarea = (Box[1]*1.0-Box[0])*(Box[3]-Box[2])
   density = n*1.0/BoxArea
   avg_dist = sqrt(1/density)
   Box = [min(x[0, *])-avg_dist, max(x[0, *])+avg_dist, $
          min(x[1, *])-avg_dist, max(x[1, *])+avg_dist]
   ;print, Box
ENDIF
IF ~keyword_set(Boundary) THEN B = ['R', 'R'] ELSE B = Boundary

IF strcmp(B[0], 'P') + strcmp(B[0], 'R') EQ 0 THEN B[0] = 'R'
IF strcmp(B[1], 'P') + strcmp(B[1], 'R') EQ 0 THEN B[1] = 'R'

Boxarea = (Box[1]*1.0-Box[0])*(Box[3]-Box[2])
density = n*1.0/BoxArea
avg_dist = sqrt(1/density)
thres = max([8*avg_dist, 8*max(w)])
IF keyword_set(threshold) THEN thres = threshold

xp = createperiodicboundary(x, Box, B)
wp = w
FOR i = 0, (kd_length(xp, 2)/n - 1) DO BEGIN
   wp = [wp, w]
ENDFOR

tmp = total(xp^2,1) - wp^2

qhull, [xp,transpose(tmp)], tr

; at least one vertex in the triangulation needs to be
; a real particle (there are 'n' real particles)
i = where(tr[0, *] LT n OR tr[1, *] LT n OR tr[2, *] LT n)
tr = tr[*, i]

; next part:  remove any triangle that contains an edge
; longer than 'thresh' --> presumably an edge going toward
; infinity
j = make_array(kd_length(tr, 2), value=0)
FOR i = 0, kd_length(tr, 2)-1 DO BEGIN
   IF max([norm(xp[*,tr[0,i]]-xp[*,tr[1,i]]),norm(xp[*,tr[0,i]]-xp[*,tr[2,i]]),norm(xp[*,tr[1,i]]-xp[*,tr[2,i]])]) LE thres THEN j[i] = 1
ENDFOR
j=where(j eq 1)
tr = tr[*, j]


v = make_array(2, kd_length(tr, 2), value=0.0)
FOR i = 0, kd_length(tr, 2)-1 DO BEGIN & $
   v[*, i]=radicalpoint(xp[*,tr[*,i]], wp[tr[*,i]])
ENDFOR


IF n_params() GE 3 THEN BEGIN

   cnt = make_array(n, value=0L)
   vorcells = 0
   
   FOR i = 0, n-1 DO BEGIN
      
      j = where(tr[0, *] EQ i OR tr[1, *] EQ i OR tr[2, *] EQ i)
      IF kd_length(j) GT 2 AND j[0] NE -1 THEN BEGIN
         qhull, v[0,j],v[1,j], order
         order2 = j
         order2[0] = order[0, 0]
         order2[1] = order[1, 0]
         FOR k = 2, kd_length(j)-1 DO BEGIN
            order2[k] = order[1, where(order[0, *] EQ order2[k-1])]
         ENDFOR
         j = j[order2]
         
         vorcells = [vorcells, j]
         cnt[i] = kd_length(j)
         IF max(tr[*, j]) GT n-1 THEN ob[i] = 1
      ENDIF
      
   ENDFOR
   
   ri = [n+1,total(cnt,/cum)+n+1,vorcells[1:kd_length(vorcells)-1]] 
   
ENDIF

tri = tr

IF arg_present(area) THEN BEGIN 
   area = make_array(kd_length(x, 2))
   FOR j = 0, kd_length(area)-1 DO BEGIN
      area[j] = poly_area(v[0, ri[ri[j]:ri[j+1]-1]], v[1, ri[ri[j]:ri[j+1]-1]])
   ENDFOR

ENDIF


return, v


END