// Demonstrates InCircle test for Delaunay Triangulation

#include <stdio.h>
#include <math.h>

#define X   0
#define Y   1
#define Z   2

#define MAX (100)

// These map coordinates to 0 ... 530
#define TRANSX(X) (530.0*((X)-xmin)/xdiff)
#define TRANSY(Y) (530.0*((Y)-ymin)/ydiff)

#define sqr(X) ((X)*(X))

/*---------------------------------------------------------------------
O'Rourke
---------------------------------------------------------------------*/
double  Volumei(int a[3], int b[3], int c[3], int p[3])
{
   double     vol;
   double  ax, ay, az, bx, by, bz, cx, cy, cz;

   ax = a[X] - p[X];
   ay = a[Y] - p[Y];
   az = a[Z] - p[Z];
   bx = b[X] - p[X];
   by = b[Y] - p[Y];
   bz = b[Z] - p[Z];
   cx = c[X] - p[X];
   cy = c[Y] - p[Y];
   cz = c[Z] - p[Z];

   vol =   ax * (by*cz - bz*cy)
         + ay * (bz*cx - bx*cz)
         + az * (bx*cy - by*cx);

   return vol;
}

int     AreaSign2( int a[3], int b[3], int c[3] )
{
    double area2;

    area2 = ( b[0] - a[0] ) * (double)( c[1] - a[1] ) -
            ( c[0] - a[0] ) * (double)( b[1] - a[1] );

    /* The area should be an integer. */
    if      ( area2 >  0.5 ) return  1;
    else if ( area2 < -0.5 ) return -1;
    else                     return  0;
}

main()
{
int pt[MAX][3],n=0;
double A,B,C,D,E,F,G;
double center[2],radius;
double xmin,ymin,xmax,ymax,xdiff,ydiff,xmid,ymid;
double vol;
int asign,i;

while (n<MAX && scanf("%d %d",&pt[n][X],&pt[n][Y]) != EOF)
{
  pt[n][Z]=sqr(pt[n][X])+sqr(pt[n][Y]);
  n++;
}
if (n==MAX)
{
  printf("Too many points %d\n",__LINE__);
  exit(0);
}
if (n<3)
{
  printf("At least 3 points are needed %d\n",__LINE__);
  exit(0);
}

// Center and radius of circle for triangle, p. 189
A=pt[1][0]-pt[0][0];
B=pt[1][1]-pt[0][1];
C=pt[2][0]-pt[0][0];
D=pt[2][1]-pt[0][1];
E=A*(pt[0][0]+pt[1][0])+B*(pt[0][1]+pt[1][1]);
F=C*(pt[0][0]+pt[2][0])+D*(pt[0][1]+pt[2][1]);
G=2.0*(A*(pt[2][1]-pt[1][1])-B*(pt[2][0]-pt[1][0]));
if (G==0.0)
{
  fprintf(stderr,"Fatal G=0 at %d\n",__LINE__);
  exit(0);
}
center[X]=(D*E-B*F)/G;
center[Y]=(A*F-C*E)/G;
radius=sqrt(sqr(pt[0][X]-center[X])+sqr(pt[0][Y]-center[Y]));

xmin=xmax=pt[0][X];
for (i=1;i<n;i++)
  if (pt[i][X]<xmin)
    xmin=pt[i][X];
  else if (pt[i][X]>xmax)
    xmax=pt[i][X];
if (center[X]+radius>xmax)
  xmax=center[X]+radius;
if (center[X]-radius<xmin)
  xmin=center[X]-radius;

ymin=ymax=pt[0][Y];
for (i=1;i<n;i++)
  if (pt[i][Y]<ymin)
    ymin=pt[i][Y];
  else if (pt[i][Y]>ymax)
    ymax=pt[i][Y];
if (center[Y]+radius>ymax)
  ymax=center[Y]+radius;
if (center[Y]-radius<ymin)
  ymin=center[Y]-radius;

xdiff=xmax-xmin;
ydiff=ymax-ymin;
//Patch to scale identically so perpendicular relationships hold
if (xdiff>ydiff)
  ydiff=xdiff;
else
  xdiff=ydiff;

printf("%%!PS\n");
printf("%%%%BoundingBox: 0 0 605 755\n");
printf("%%%%EndComments\n");
printf("40 155 translate\n");
printf("newpath\n");
printf("%f %f 8 0 360 arc stroke\n", TRANSX(center[X]),TRANSY(center[Y]));
printf("%f %f %f 0 360 arc stroke\n", TRANSX(center[X]),TRANSY(center[Y]),
       TRANSX(center[X]+radius)-TRANSX(center[X]));

asign=AreaSign2(pt[0],pt[1],pt[2]); // Just in case not counterclockwise
fprintf(stderr,"asign is %d\n",asign);

for (i=0;i<n;i++)
{
  vol=Volumei(pt[0],pt[1],pt[2],pt[i]);
  fprintf(stderr,"%d %d gives vol of %f\n",pt[i][X],pt[i][Y],vol);
  if (vol*asign==0.0) // On circle
    printf("%f %f 4 0 360 arc stroke\n",TRANSX(pt[i][X]),TRANSY(pt[i][Y]));
  else if (vol*asign>0.0) // Inside circle
    printf("%f %f 2 0 360 arc fill\n",TRANSX(pt[i][X]),TRANSY(pt[i][Y]));
  else // Outside circle
    printf("%f %f 6 0 360 arc fill\n",TRANSX(pt[i][X]),TRANSY(pt[i][Y]));
}
printf("showpage\n\n");
printf("%%EOF\n");
}