/****************************************************************************
 * Implementation of Preparata and Hong's divide-and-conquer algorithm      *
 * for the computation of convex hulls in two dimension                     *
 * Written by: Abhijit Das (abhij@cse.iitkgp.ernet.in)                      *
 * Last modified on Aug 14, 2006                                            *
 ****************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

/* default size of the data set, overridden by command line argument */
#define N_DEFAULT 100

/* data type: point (in two-dimension) */
typedef struct {
   double x;        /* x-coordinate */
   double y;        /* y-coordinate */
} point;

/* data type: convex hull */
typedef struct {
   int size;        /* number of vertices on the hull */
   point *vertex;   /* list of vertices on the hull in the clockwise order
                       starting from the leftmost vertex */
} hull;

/* Generate a list of n random points */
point *genInput ( int n )
{
   point *P = NULL;
   int i;

   /* seed the random number generator by current time */
   srand((unsigned int)time(NULL));

   P = (point *)malloc(n * sizeof(point));
   for (i=0; i<n; ++i) {
      P[i].x = (double)rand() / 2147483648.;
      P[i].y = (double)rand() / 2147483648.;
   }
   return P;
}

/* merge-sort a list of points in increasing order of x-coordinate */
void sortInput ( point *P, int n )
{
   int i,j,k,n1,n2;
   point *Q;

   if (n == 1) return;
   n1 = n/2; n2 = n-n1;
   sortInput(P,n1);
   sortInput(&P[n1],n2);
   Q = (point *)malloc(n*sizeof(point));
   i = 0; j = n1;
   for (k=0; k<n; ++k) {
      if (i>=n1) Q[k] = P[j++];
      else if (j>=n) Q[k] = P[i++];
      else if (P[i].x < P[j].x) Q[k] = P[i++];
      else Q[k] = P[j++];
   }
   for (k=0; k<n; ++k) P[k] = Q[k];
   free(Q);
}

/* print the input points to the data file "CH1.dat" */
void printInput ( point *P, int n )
{
   int i;
   FILE *fp;

   fp = (FILE *)fopen("CH1.dat","w");
   for (i=0; i<n; ++i) {
      fprintf(fp,"%12.10lf %12.10lf\n", P[i].x, P[i].y);
   }
   fclose(fp);
}

/* returns true if and only if P3 lies to the left of the oriented line
   segment P1P2 */
int toLeft ( point P1, point P2, point P3 )
{
   return( (P1.x * P2.y - P1.y * P2.x) +
           (P2.x * P3.y - P2.y * P3.x) +
           (P3.x * P1.y - P3.y * P1.x)
           > 0 );
}

/* returns true if and only if P3 lies to the right of the oriented line
   segment P1P2 */
int toRight ( point P1, point P2, point P3 )
{
   return( (P1.x * P2.y - P1.y * P2.x) +
           (P2.x * P3.y - P2.y * P3.x) +
           (P3.x * P1.y - P3.y * P1.x)
           < 0 );
}

/* returns true if and only if P1P2 is an upper tangent to the hull H
   touching at the j-th vertex */
int upperTangent ( hull H, int j, point P1, point P2 )
{
   int i,k;

   i = (j==0) ? H.size-1 : j-1;  /* the previous vertex */
   k = (j==H.size-1) ? 0 : j+1;  /* the next vertex */
   return (toRight(P1,P2,H.vertex[i]) && toRight(P1,P2,H.vertex[k]));
}

/* returns true if and only if P1P2 is a lower tangent to the hull H
   touching at the j-th vertex */
int lowerTangent ( hull H, int j, point P1, point P2 )
{
   int i,k;

   i = (j==0) ? H.size-1 : j-1;  /* the previous vertex */
   k = (j==H.size-1) ? 0 : j+1;  /* the next vertex */
   return (toLeft(P1,P2,H.vertex[i]) && toLeft(P1,P2,H.vertex[k]));
}

/* merge the left and the right hulls */
hull mergeHull ( hull LH, hull RH )
{
   int i, j, l1, l2, r1, r2;
   hull H;

   /* l1 and l2 should point to the rightmost vertex of LH, and
      r1 and r2 should point to the leftmost vertex of RH */
   l1 = 0;
   while ((l1 <= LH.size-2) && (LH.vertex[l1+1].x > LH.vertex[l1].x)) ++l1;
   l2 = l1;
   r1 = r2 = 0;

   /* compute the common upper tangent */
   while ( (!upperTangent(LH,l1,LH.vertex[l1],RH.vertex[r1])) ||
           (!upperTangent(RH,r1,LH.vertex[l1],RH.vertex[r1]))
         ) {
      while (!upperTangent(LH,l1,LH.vertex[l1],RH.vertex[r1])) --l1;
      while (!upperTangent(RH,r1,LH.vertex[l1],RH.vertex[r1])) ++r1;
   }

   /* compute the common lower tangent */
   while ( (!lowerTangent(LH,l2,LH.vertex[l2],RH.vertex[r2])) ||
           (!lowerTangent(RH,r2,LH.vertex[l2],RH.vertex[r2]))
         ) {
      while (!lowerTangent(LH,l2,LH.vertex[l2],RH.vertex[r2])) {
         l2 = (l2 == LH.size-1) ? 0 : l2+1;
      }
      while (!lowerTangent(RH,r2,LH.vertex[l2],RH.vertex[r2])) {
         r2 = (r2 == 0) ? RH.size-1 : r2-1;
      }
   }

   /* the maximum size of the merge hull is the sum of the sizes of LH and RH */
   H.vertex = (point *)malloc((LH.size + RH.size) * sizeof(point));

   j = 0; /* j is used for writing to the vertex list of the merged hull */

   /* first copy vertices from LH between indices 0 and l1 */
   i = 0;
   while (i <= l1) H.vertex[j++] = LH.vertex[i++];

   /* then copy vertices from RH between indices r1 and r2 */
   i = r1;
   while (i != r2)  {
      H.vertex[j++] = RH.vertex[i++];
      if (i == RH.size) i = 0;
   }
   H.vertex[j++] = RH.vertex[r2];

   /* finally copy vertices from LH between indices l2 and LH.size-1 */
   i = l2;
   while (i != 0) {
      H.vertex[j++] = LH.vertex[i++];
      if (i == LH.size) i = 0;
   }

   /* now j stores the actual size of the merged hull */
   H.size = j;

   return H;
}

/* Preparata and Hong's divide-and-conquer algorithm */
hull genHull ( point *P, int n )
{
   hull H, LH, RH;
   int n1, n2;

   if (n == 1) {         /* convex hull on one vertex */
      H.size = 1;
      H.vertex = (point *)malloc(sizeof(point));
      H.vertex[0] = P[0];
   } else if (n == 2) {  /* convex hull on two vertices */
      H.size = 2;
      H.vertex = (point *)malloc(2 * sizeof(point));
      H.vertex[0] = P[0];
      H.vertex[1] = P[1];
   } else if (n == 3) {  /* convex hull on three vertices */
      H.size = 3;
      H.vertex = (point *)malloc(3 * sizeof(point));
      H.vertex[0] = P[0];
      if (toRight(P[0],P[1],P[2])) {
         H.vertex[1] = P[1];
         H.vertex[2] = P[2];
      } else {
         H.vertex[1] = P[2];
         H.vertex[2] = P[1];
      }
   } else {              /* convex hull on four or more vertices */
      n1 = n/2;
      n2 = n - n1;
      LH = genHull(P,n1);       /* recursively generate the left hull */
      RH = genHull(&P[n1],n2);  /* recursively generate the right hull */
      H = mergeHull(LH,RH);     /* merge the two sub-hulls */
      free(LH.vertex);
      free(RH.vertex);
   }
   return H;
}

/* print the computed convex hull to the data file "CH2.dat" */
void printOutput ( hull H )
{
   int i;
   FILE *fp;

   fp = (FILE *)fopen("CH2.dat","w");
   for (i=0; i<H.size; ++i) {
      fprintf(fp,"%12.10lf %12.10lf\n", H.vertex[i].x, H.vertex[i].y);
   }
   fprintf(fp,"%12.10lf %12.10lf\n", H.vertex[0].x, H.vertex[0].y);
   fclose(fp);
}

int main ( int argc, char *argv[] )
{
   int n;
   point *P;
   hull H;

   /* argv[1] overrides the default size of the input point set */
   n = (argc > 1) ? atoi(argv[1]) : N_DEFAULT;

   P = genInput(n);      /* generate the input point set */
   printf("%d random points generated...\n", n);
   sortInput(P,n);       /* sort the input point set w.r.t. x-coordinates */
   printf("Input points sorted...\n");
   printInput(P,n);      /* print input points to data file "CH1.dat" */
   H = genHull(P,n);     /* generate the convex hull recursively */
   printf("Convex hull generated: size = %d...\n", H.size);
   printOutput(H);       /* print output hull to data file "CH2.dat" */
   free(P);
}