; $Id: voronoi.pro,v 1.2 1993/08/20 14:40:44 dave Exp $

function isright, x0, y0, x1, y1, x2, y2 
; return 1 if Pnt0 is to right of Pnt1-> Pnt2
; return 0 if it is on the line.
; return -1 if Pnt0 is to the left of Pnt1 -> Pnt2

z = (x0-x1) * (y2-y1) - (y0-y1) * (x2-x1)
if z gt 0.0 then return, 1
if z lt 0.0 then return, -1
return, 0
end


PRO voronoi, x, y, i0, c, xp, yp
; Copyright (c) 1992, Research Systems, Inc.  All rights reserved.
;	Unauthorized reproduction prohibited.
;+
; NAME:
;	VORONOI
; PURPOSE:
;	Compute the Voronoi polygon of a point within an irregular
;	grid of points, given the Delaunay triangulation.  The Voronoi
;	polygon of a point contains the region closer to that point
;	than to any other point.
; CATEGORY:
;	Gridding.
; CALLING SEQUENCE:
;	VORONOI, x, y, i0, c, xp, yp
; INPUTS:
;	x = array containing X locations of the points.
;	y = y array of point coordinates.
;	i0 = index of point.
;	c = connectivity list from the Delaunay triangulation.
;		This list is produced with the CONNECTIVITY keyword
;		of the TRIANGULATE procedure.
; OUTPUTS:
;	xp, yp = vertices of voroni polygon.
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;	None.
; RESTRICTIONS:
;	The polygons only cover the convex hull of the set of points.
; PROCEDURE:
;	For interior points, the polygon is constructed by connecting
;	the midpoints of the lines connecting the point with its Delaunay
;	neighbors.  Polygons are traversed in a counterclockwise direction.
;
;	For exterior points, the set described by the midpoints of the
;	connecting lines, plus the circumcenters of the two triangles
;	that connect the point to the two adjacent exterior points.
; EXAMPLE:
;	To draw the voroni polygons of each point of an irregular 
;	grid:
;	x = randomu(seed, n)		;Random grid of N points
;	y = randomu(seed, n)
;	triangulate, x, y, tr, CONN=c	;Triangulate it
;	for i=0, n-1 do begin
;		voronoi, x, y, i, c, xp, yp	;Get the ith polygon
;		polyfill, xp, yp, color = (i mod 10) + 2  ;Draw it
;		endfor
;
; MODIFICATION HISTORY:
;	DMS, RSI.	Dec, 1992.  Original version.
;-


p = c(c(i0):c(i0+1)-1)  	;Verts of polygon
np = n_elements(p)

if i0 ne p(0) then begin	;Boundary point?
    xp = fltarr(np, /NOZERO)	;No. interior point.
    yp = fltarr(np, /NOZERO)
    for i=0, np-1 do begin
	m = p(i)
	j = p((i + 1) mod np)	;Successor
	cir_3pnt, x([m,j,i0]), y([m,j,i0]), r, x0, y0
	xp(i) = x0
	yp(i) = y0
	endfor
endif else begin		;Boundary point..
    q = [p(0), p(1), p(2)]
    cir_3pnt, x(q), y(q), r, x0, y0 ;first triangle
    x1 = (x(p(1)) + x(i0))/2.	;First midpoint
    y1 = (y(p(1)) + y(i0))/2.
    l1 =  isright(x1, y1, x(i0), y(i0), x0, y0) eq 1	;True to add midpoint
    q = [p(0), p(np-1), p(np-2)]	;Last triangle
    cir_3pnt, x(q), y(q), r, x2, y2
    x3 = (x(p(np-1)) + x(i0))/2.	;Last midpoint
    y3 = (y(p(np-1)) + y(i0))/2.
    l2 = isright(x3, y3, x(i0), y(i0), x2, y2) eq -1  ;True to add last midpnt
    n = np-1 + l1 + l2		;# of vertices
    xp = fltarr(n, /NOZERO)
    yp = fltarr(n, /NOZERO)
    xp(0) = x0			;Add first triangle's circumcenter
    yp(0) = y0
    for i=2, np-3 do begin	;Interior triangles
	q = [ i0, p(i), p(i+1) ]
	cir_3pnt, x(q), y(q), r, x0, y0
	xp(i-1) = x0
	yp(i-1) = y0
	endfor
    xp(np-3) = x2		;Last triangle's circumcenter
    yp(np-3) = y2
    l = np-2
    if l2 then begin		;Add last midpnt?
	xp(l) = x3
	yp(l) = y3
	l = l + 1
	endif
    xp(l) = x(i0)		;Add central point
    yp(l) = y(i0)
    l = l + 1
    if l1 then begin		;Add first midpnt?
	xp(l) = x1
	yp(l) = y1
	l = l + 1
	endif
    endelse
end