// chull_graham.H
// find the convex hull of a set of 2d points (C++ template)
// Graham scan: Sort points and add to convex hull in sorted order
// algorithm by Andrew (IPL, 1979) based on algorithm by Graham (IPL, 1972)
// Version : 1.0
// Last modified: 1999/08/23 18:49:00
// Author : R. Wenger

/* Copyright (c) 1998 by The Ohio State University
 *
 * Permission to use, copy, modify, distribute, and sell this software 
 * and its documentation for any purpose is hereby granted without fee,
 * provided that the above copyright notice appear in all copies and/or
 * modified versions and that both that copyright notice and this 
 * permission notice appear in supporting documentation.  This software 
 * is provided "as is" without express or implied warranty, including 
 * but not limited to the warranties of merchantibility, fitness for 
 * a particular purpose and noninfringement.
 *
 */

#ifndef CHULL_GRAHAM_H
#define CHULL_GRAHAM_H

#include <algo.h>


#ifndef SWAP_ARG
#define SWAP_ARG

// function adaptor. swap function arguments.
// swap_arg(func)(a,b) = func(b, a)

template <class ARG>
  class swap_arg {

    protected :
      bool (*F)(ARG &, ARG &);

    public :
      swap_arg(bool (*func) (ARG &, ARG &)) : 
	F(func) {};
      bool operator()(ARG & a, ARG & b) const 
        { return(F(b,a)); }
  };

#endif SWAP_ARG


template <class POINT_ITER, class POINT, class ORESULT>
void chull_graham(POINT_ITER firstp, POINT_ITER endp, 
		  bool (*compare) (POINT &, POINT &), 
		  ORESULT (*orient) (POINT &, POINT &, POINT &), 
		  POINT_ITER & endc)
// reorder point array [firstp,endp) so that the points [firstp,endc) are 
// distinct, extremal vertices of the convex hull in counter-clockwise order
//
// POINT_ITER = type of random access iterator to point array
// POINT = point type. POINT_ITER = (POINT *)
// ORESULT = type returned by function orient
// firstp = iterator (pointer) to first element in array [firstp,endp)
// endp = iterator (pointer) to past-the-end of array [firstp,endp)
// orient = calculate orientation
// compare = point comparison function
// endc = iterator (pointer) to past-the-end of the convex hull points

/* 
   compare(p1,p2) = 
     return 1 if p1.x < p2.x or (p1.x == p2.x and p1.y < p2.y)
     return 0 otherwise
   orient(p1,p2,p3) = returns orientation of points p1, p2, p3 
     (pos = left, neg = right, 0 = collinear)
*/

{
  const POINT_ITER lastp = endp - 1;  // iterator to last point

  if (endp - firstp < 2) 
    { endc = endp; return; };

// at least two points in [firstp, endp)

/* find min/max points and put in first/last position */
  POINT_ITER imin = min_element(firstp, endp, compare);
  if (imin != firstp)
    iter_swap(imin, firstp);

  POINT_ITER imax = max_element(firstp, endp, compare);

  if (!compare(*firstp, *imax))
/* all points have same x/y coordinates */
    { endc = firstp+1; return; };  

  iter_swap(imax, lastp);

// divide points into those above and those below segment from
//   point *firstp to point *lastp
  POINT_ITER liter = firstp + 1;  /* left iterator */
  POINT_ITER riter = lastp - 1;   /* right iterator */

  while (liter <= riter) {
    if (orient(*firstp,*lastp,*liter) <= 0)
      liter++;
    else if (orient(*firstp,*lastp,*riter) > 0)
      riter--;
    else
      { iter_swap(liter,riter); liter++; riter--; };
    };
  
  iter_swap(liter,lastp);
  imax = liter;

/* find the lower convex hull */

  POINT_ITER lastc = firstp + 1;

/* sort first half of the table in increasing order */
  if (lastc < imax)
    sort(lastc, imax, compare);

  POINT_ITER curp = firstp + 2;

  while (curp <= imax) {
    while (lastc != firstp && orient(*(lastc-1), *lastc, *curp) <= 0)
      lastc--;

/* add *curp to list of vertices */
    lastc++;
    if (lastc < curp)
      iter_swap(lastc, curp);

    curp++;
  };

  imax = lastc;

/* find the upper convex hull */

/* sort second half of the table in decreasing order */
  if (curp < lastp)
    sort(curp, endp, swap_arg<POINT>(compare));

  while (curp < endp) {
    while (lastc != imax && orient(*(lastc-1), *lastc, *curp) <= 0)
      lastc--;

/* add *curp to list of vertices */
    lastc++;
    if (lastc < curp)
      iter_swap(lastc, curp);

    curp++;
  };

/* eliminate any trailing non-extreme points */
  while (lastc != imax && orient(*(lastc-1), *lastc, *firstp) <= 0)
    lastc--;

  endc = lastc+1;
};

#endif CHULL_GRAHAM_H