Viewing contents of file '../idllib/jhuapls1r/usr/convexhull.pro'
;-------------------------------------------------------------
;+
; NAME:
;       CONVEXHULL
; PURPOSE:
;       Return the convex hull of a polygon.
; CATEGORY:
; CALLING SEQUENCE:
;       convexhull, x, y, xh, yh
; INPUTS:
;       x,y = original polygon vertices.       in
; KEYWORD PARAMETERS:
; OUTPUTS:
;       xh,yh = convex hull polygon vertices.  out
; COMMON BLOCKS:
; NOTES:
;       Notes: The convex hull of a polygon is the minimum polygon
;         that circumscribes the original polygon.  It is the shape
;         a rubber band would take if placed around the original
;         polygon.
; MODIFICATION HISTORY:
;       R. Sterner, 2 Oct, 1990
;       R. Sterner, 26 Feb, 1991 --- renamed from convex_hull.pro
;
; Copyright (C) 1990, Johns Hopkins University/Applied Physics Laboratory
; This software may be used, copied, or redistributed as long as it is not
; sold and this copyright notice is reproduced on each copy made.  This
; routine is provided as is without any express or implied warranties
; whatsoever.  Other limitations apply as described in the file disclaimer.txt.
;-
;-------------------------------------------------------------
 
	function onhull, x1, y1, x2, y2, x, y
 
	;--- (x1,y1) and (x2,y2) are the two endpoints of a polygon side.
	;--- x,y are the arrays of the polygon vertices.
	ux = y2 - y1		  ; This is a vector perpendicular to
	uy = x1 - x2		  ; the polygon side.
	d = ux*(x-x2) + uy*(y-y2) ; Dot prod of U vect to all vertices.
	s = sign(d)		  ; Side of polygon side each vertex is on.
	r = max(s) - min(s)	  ; 1: one side, -1: other side, 0: bewteen.
	return, r eq 1		  ; Seg on C.H. if all pts on one side.
	end
 
 
;-------  Convexhull.pro = return the convex hull of a polygon  ------
;	R. Sterner, 3 Oct, 1990
 
	pro convexhull, x, y, xh, yh, help=hlp
 
	if (n_params(0) lt 4) or keyword_set(hlp) then begin
	  print,' Return the convex hull of a polygon.'
	  print,' convexhull, x, y, xh, yh'
	  print,'   x,y = original polygon vertices.       in'
	  print,'   xh,yh = convex hull polygon vertices.  out'
	  print,' Notes: The convex hull of a polygon is the minimum polygon'
	  print,'   that circumscribes the original polygon.  It is the shape'
	  print,'   a rubber band would take if placed around the original'
	  print,'   polygon.'
	  return
	endif
 
;---  Conjecture: The mean of the polygon points is inside the convex hull ----
	xm = mean(x)
	ym = mean(y)
;---  Conjecture: The farthest point from any point inside the convex hull ----
	;                 is on the convex hull.
	d = (x-xm)^2 + (y-ym)^2		; Dist^2 from mean to all polygon pts.
	w = where(d eq max(d))		; Find farthest point.
	i0 = w(0)			; Index of farthest point.
	xs = shift(x,-i0)		; Shift arrays to put farthst pt first.
	ys = shift(y,-i0)
	n = n_elements(x)		; Size of polygon.
	xh = fltarr(1) + xs(0)		; Start convex hull array.
	yh = fltarr(1) + ys(0)
	l = 0				; Last convex hull point index.
 
	i0 = 1				; Search start index.
 
	;-----  Starting with a known convex hull point find the next  -------
loop:
	for i = i0, n-1 do begin
	  if onhull(xh(l),yh(l),xs(i),ys(i),xs,ys) then begin
	    xh = [xh,xs(i)]		; Add new point to hull.
	    yh = [yh,ys(i)]
	    l = l + 1			; Count hull point.
	    i0 = i + 1			; Start search at next point.
	    goto, loop
	  endif
	endfor
 
	return
 
	end