/*****
** Duesseldorf, den 10.4.1995.
** chull.C: Berechnung der Konvexen Huelle von Punkten im R^2.
** Wenn die Eingabedaten exakt in einfacher Genauigkeit darstellbar sind, dann 
** arbeitet chull ohne Rundungsfehler.
*****/

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

#include "fpu_ctrl.h"
#include "sgnsum.h"	// ESSA
#include "interval.h"	// Intervall-Paket mit einfacher Genauigkeit.

/*----------------------------------------------------------------------------
type tStat:
Beschreibung: Struktur zur Aufnahme der Statistik.
----------------------------------------------------------------------------*/
typedef struct
{
 int n,			// Anzahl der Punkte
     nHull,		// Anzahl der Punkte der berechneten Huelle.
     nGraham,		// Anzahl der nach dem Loeschen an Graham ueberg. Pkt.
     time,		// Gesamtlaufzeit.
     nSgnsumPre,	// Anzahl der Aufrufe von ESSA im Preprocessing.
     nSgnsumGR,		// Anzahl der Aufrufe von ESSA in Graham-Scan.
     timeSgnsumPre,	// Laufzeit fuer ESSA im Preprocessing.
     timeSgnsumGR;	// Laufzeit fuer ESSA in Graham-Scan.
} tStat;

tStat st; // globale Statistikvariable.

/*----------------------------------------------------------------------------
class cPoint:
Beschreibung: Darstellung der Punkte in R^2.
----------------------------------------------------------------------------*/
class cPoint
{
 public:
  cPoint();
  cPoint(float,float);
  cPoint(cPoint&);
  cPoint &operator=(cPoint&);
 
  float x,y;
};

cPoint::cPoint()
{
 x=y=0;
}

cPoint::cPoint(float _x, float _y)
{
 x= _x;
 y= _y;
}

cPoint::cPoint(cPoint &p)
{
 x= p.x;
 y= p.y;
}

cPoint &cPoint::operator=(cPoint &p)
{
 x= p.x;
 y= p.y;
 return *this;
}

/*----------------------------------------------------------------------------
class cIntStack:
Beschreibung: Realisierung eines Stacks zur Aufnahme von Integern.
----------------------------------------------------------------------------*/
typedef struct sIntStack
{
 int i;
 struct sIntStack *next;
}tIntStack;

class cIntStack
{
 public:
  cIntStack();		// Initialisierung.
  ~cIntStack();		// Destructor.
  int IsEmpty(void);	// Test, ob Stack leer ist.
  void Push(int p);	// Einfuegen in Stack.
  void Pop(void);	// Erstes Element aus Stack entfernen.
  int Get(int n=0);	// n=0: Erstes Element lesen, n=1 zweites Element lesen.
 
  tIntStack *top;
};

cIntStack::cIntStack()
{
 top=NULL;
}


cIntStack::~cIntStack()
{
 while (!IsEmpty()) Pop();
}

cIntStack::IsEmpty(void)
{
 return top==NULL;
}

void cIntStack::Push(int p)
{
 tIntStack *IS;
 
 IS= new tIntStack;
 IS->i= p;
 IS->next= top;
 top= IS;
}

void cIntStack::Pop(void)
{
 tIntStack *IS;
 
 if (top!=NULL)
 {
  IS=top;
  top= top->next;
  free(IS);
 }
}

int cIntStack::Get(int n)
{
 if (n==0) return top->i;
 else return top->next->i;
}

/*----------------------------------------------------------------------------
Funktion swap:
Beschreibung: Vertauscht a und b.
----------------------------------------------------------------------------*/
void swap(cPoint &a, cPoint &b)
{
 cPoint tmp= a;
 a= b;
 b= tmp;
}

/*----------------------------------------------------------------------------
Funktion LeftTest:
Input: a,b,c: cPoint, ESSA..Anzahl der ESSA-Aufrufe, tESSA..Laufzeit der 
   ESSA-Aufrufe
Output: 1 oder 0.
Beschreibung: Es wird getestet, ob die Drehung von b nach y im Ursprung 
   a eine Linksdrehung ist. Fuehrt der Test mit Hilfe von Intervallarithmetik 
   nicht zum Ergebnis, wird ESSA aufgerufen. Dies garantiert immer eine
   exakte Lösung.
----------------------------------------------------------------------------*/
int LeftTest(cPoint a, cPoint b, cPoint y, int &ESSA, int &tESSA)
{
 interval d,Si[9];
 int i,CPUT;
 double S[9];

//** Die Summanden: 
 Si[1]= interval(y.x)*interval(a.y);
 Si[2]= interval(a.x)*interval(b.y);
 Si[3]= interval(b.x)*interval(y.y);
 Si[4]= -interval(y.x)*interval(b.y);
 Si[5]= -interval(a.x)*interval(y.y);
 Si[6]= -interval(b.x)*interval(a.y);
 
 d=0;
 for (i=1;i<=6;i++) d+= Si[i];
      
// printf("d: %e %e, w(d): %e\n",d.Inf(),d.Sup(),d.Diam());
 
 if (d.Inf()>0) return 1;
 if (d.Sup()<0) return -1;
 if (d.Inf==0&&d.Sup==0) return 0;

//** Vorbereitungen fuer den Aufruf von ESSA. 
 S[1]= (double)y.x*(double)a.y;
 S[2]= (double)a.x*(double)b.y;
 S[3]= (double)b.x*(double)y.y;
 S[4]= -(double)y.x*(double)b.y;
 S[5]= -(double)a.x*(double)y.y;
 S[6]= -(double)b.x*(double)a.y;
  
 ESSA++;
 CPUT=-GetCPUUsage();
// printf("sgnsum %i\n",st.nSgnsum);
 i=sgnsum(S,6);
 CPUT+=GetCPUUsage();
 tESSA+=CPUT;
 return i;
}

/*----------------------------------------------------------------------------
Funktion FindLowest:
Beschreibung: Suche nach dem am weitesten unten rechts gelegenen Punkt in P.
   Dieser wird nach P[0] kopiert.
----------------------------------------------------------------------------*/
//** FindLowest.
cPoint FindLowest(int &n, cPoint *P)
{
 int i,m,j;
 
 m=0;
 
 for (i=1;i<n;i++)
  if ( (P[i].y<P[m].y) || 
       (P[i].y==P[m].y && P[i].x>P[m].x) )
   m=i;
 swap(P[m],P[0]);
 j=1;
 for (i=1;i<n;i++) if (P[i].x!=P[0].x||P[i].y!=P[0].y) P[j++]=P[i];
 n= j;
}

/*----------------------------------------------------------------------------
Funktion SortCounterClockwise:
Input: n..Laenge der Liste ra, ra: Liste von cPoint, P0: Der Ursprung.
Beschreibung: Die Elemente aus ra werden entgegen dem Uhrzeigersinn mit 
   Ursprung P0 sortiert. Der Algorithmus basiert auf Heapsort.
----------------------------------------------------------------------------*/
void SortCounterClockwise(int n, cPoint *ra, cPoint P0)
{
 int l,j,ir,i;
 cPoint rra;
 
 l=(n >> 1)+1;
 ir=n;
 for(;;)
 {
  if (l>1)
   rra=ra[--l];
  else
  {
   rra=ra[ir];
   ra[ir]=ra[1];
   if (--ir==1)
   {
    ra[1]=rra;
    return;
   }
  } /* (l<=1) */
  i=l;
  j=l<<1;
  while (j<=ir)
  {
   if ( j<ir && LeftTest(P0,ra[j],ra[j+1],st.nSgnsumPre,st.timeSgnsumPre)>0 ) ++j;
   if ( LeftTest(P0,rra,ra[j],st.nSgnsumPre,st.timeSgnsumPre)>0 )
   {
    ra[i]=ra[j];
    j+=(i=j);
   }
   else
    j=ir+1;
  }
  ra[i]=rra;
 } /* for(;;) */
}

/*----------------------------------------------------------------------------
Funktion DeleteDoubles:
Beschreibung: Wenn zwei Punkte auf einer Geraden durch P[0] liegen, dann wird
   der jeweils P[0] naehere Punkt geloescht.
----------------------------------------------------------------------------*/
void DeleteDoubles(int &n, cPoint *P)
{
 int i,j=1,k,CPUT;
 interval Si[13],d;
 cPoint a,b,c;
 
// printf("nDelete: %i\n",n);
 a=P[0];
 for (i=2;i<n;i++)
 {
  b= P[j];
  c= P[i];
//  printf("i: %i, j: %i\n",i,j);
  if ((b.x==c.x) && (b.y==c.y)) continue;
  if (LeftTest(a,b,c,st.nSgnsumPre,st.timeSgnsumPre)==0)
  {
//   printf("iP[%i]: %e %e\n",i,c.x,c.y);
//   printf("jP[%i]: %e %e\n",j,b.x,b.y);
   if (b.y==c.y)
   {
    if (b.x>c.x) P[j]= P[i];
   }
   else
    if (b.y<c.y) P[j]= P[i];
  }
  else
   P[++j]= P[i];
 }
 n=j+1;
// printf("nGraham: %i\n",n);
}

/*----------------------------------------------------------------------------
Funktion Graham:
Beschreibung: Berechnung der konvexen Huelle der Punkte in P (Laenge n) mit 
   dem Graham-Scan-Algorithmus und verbessertem Leftturntest. Das Ergebnis
   wird in dem Stack ST zurückgegeben.
----------------------------------------------------------------------------*/
void Graham( cIntStack &ST, cPoint *P, int n)
{
 int i;

//** Initialize stack to P[n-1], P[0]. 
 ST.Push(n-1);
 ST.Push(0);
 if (n<=2) return;
 
//** Bottom two elements will never be removed.
 i=1;
 while (i<n)
 {
//  printf("Top: %i %i\n",ST.Get(1),ST.Get(0));
  if ( LeftTest(P[ST.Get(1)],P[ST.Get(0)],P[i],st.nSgnsumGR,st.timeSgnsumGR)>0 )
  {
   ST.Push(i);
   i++;
  }
  else ST.Pop();
 }
}

/*----------------------------------------------------------------------------
Funktion ComputeHull:
Beschreibung: Berechnung der konvexen Huelle der Punkte in der Liste P mit
   der Laenge n. Die Liste Hull der Laenge m enthaelt die entgegen dem
   Uhrzeigersinn sortierten Punkte der Huelle.
----------------------------------------------------------------------------*/
void ComputeHull(cPoint *P, int n, cPoint *Hull, int &m)
{
 cIntStack ST;
 int i;
 
//** Initialisierung der Statistik.
 st.nSgnsumPre=0;
 st.nSgnsumGR=0;
 st.nHull= 0;
 st.timeSgnsumPre=0;
 st.timeSgnsumGR=0;
 st.time=-GetCPUUsage();
 st.n= n;

//** Suche der am weitesten rechts unten gelegenen Punkt.
 FindLowest(n,P);
// printf("Lowest: %e, %e\n",P[0].x, P[0].y);

//** Sortiere die Punkte gegen den Uhrzeigersinn mit Ursprung P[0].
 SortCounterClockwise(n-1,P,P[0]);
// printf("Sorted List:\n");
// for (i=0;i<n;i++) printf("P[%i]: (%e, %e)\n",i,P[i].x,P[i].y);

//** Loesche doppelte.
 DeleteDoubles(n,P);
//printf("Delete doubles:\n");
// for (i=0;i<n;i++) printf("P[%i]: (%e, %e)\n",i,P[i].x,P[i].y);
 printf("Graham\n");
 st.nGraham= n;

//** Nun kann Graham-Scan aufgerufen werden, um die konvexe Huelle zu
//** berechnen.
 Graham(ST,P,n);

//** Das letzte Element ist doppelt. Es wird aber in der Liste gelassen, 
//** damit gnuplot ein vernuenftiges Bild liefert.
// ST.Pop();

//** Ermittlung der Gesamtlaufzeit.
 st.time+= GetCPUUsage();

 printf("GrahamEnde\n");

//** Einfuegen der Punkte aus dem Stack ST in das Array Hull.
 m=0;
 while (!ST.IsEmpty()) 
 {
//  printf("Hull[%i]: %i\n",m,ST.Get());
  Hull[m++]= P[ST.Get()];
  ST.Pop();
 } 

//** In der Statistik erscheint die richtige Anzahl von Punkten in der
//** Huelle.
 st.nHull= m-1;
}

/*----------------------------------------------------------------------------
Funktion main:
Beschreibung: Aufruf mit "chull <Dateiname>". An Dateiname wird die 
   Endung ".points" angehaengt. In der Datei "Dateiname.points" sollten
   die Punkte zu dennen die Huelle berechnet werden soll stehen.
   Dabei muss in jeder Zeile ein Punkt stehen. Die Huelle wird in der Datei
   "Dateiname.hull" gespeichert. Die Statistik wird in der Datei
   "Dateiname.stat" gespeichert.
----------------------------------------------------------------------------*/
int main(int al, char **a)
{
 cPoint P[10010],Hull[10010];
 int i,n,m;
 char c[1000],cc[1000],*ccc;
 cPoint P0,P1,P2;
 float x,y;
 FILE *f;

 if (al>2)
 {
  switch(**(a+1))
  {
   case 't':
//** Einlesen der Punkte im Textformat
    strcpy(c,*(a+2));
    strcat(c,".points");
    if ( (f=fopen(c,"rt"))==NULL )
    {
     printf("Die Datei %s konnte nicht geoeffnet werden.\n",c);
     exit(1);
    }
    n=0;
    for (;;)
    {
     if (fgets(c,900,f)==NULL) break;
     ccc= strstr(c," ");
     ccc++;
     P[n]= cPoint((float)atof(c),(float)atof(ccc));
     n++;
    }
    fclose(f);
    printf("%i Punkte gelesen\n",n);
    break;
   case 'b':

//** Einlesen der Punkte im Binaerformat
    strcpy(c,*(a+2));
    strcat(c,".bpoints");
    if ( (f=fopen(c,"r"))==NULL )
    {
     printf("Die Datei %s konnte nicht geoeffnet werden.\n",c);
     exit(1);
    }
    fread(&n,1,sizeof(int),f);
    for (i=0;i<n;i++)
    {
     if (fread(&x,1,sizeof(float),f)==0) break;
     if (fread(&y,1,sizeof(float),f)==0) break;
     P[i].x= x;
     P[i].y= y;
    }
    fclose(f);
    printf("%i Punkte gelesen\n",n);
//** Speicherung der Punkte im Textformat fuer gnuplot.
    RmDefault();
    strcpy(c,*(a+2));
    strcat(c,".points");
    f= fopen(c,"wt");
    for (i=0;i<n;i++) fprintf(f,"%1.8e %1.8e\n",P[i].x,P[i].y);
    fclose(f);
    break;
   default:
    printf("b fuer binaere Dateien, t fuer Textdateien\n");
  }

//** Berechnung der Huelle.
    printf("Starte ComputeHull\n"); 
    ComputeHull(P,n,Hull,m);

//** Speicherung der Huelle
    RmDefault();
    strcpy(c,*(a+2));
    strcat(c,".hull");
    f= fopen(c,"wt");
    for (i=0;i<m;i++) fprintf(f,"%1.8e %1.8e\n",Hull[i].x,Hull[i].y);
    fclose(f);


//** Speicherung der Statistik.
    strcpy(c,*(a+2));
    strcat(c,".stat");
    f= fopen(c,"wt");
    fprintf(f,"%5i & %5i & %5i & %5i & %5i & %5i & %5i & %5i\\\\\n",
            st.n,st.nHull,st.nGraham,st.nSgnsumPre,st.timeSgnsumPre,
            st.nSgnsumGR,st.timeSgnsumGR,st.time);
    fprintf(f,"  n   | nHull |  nGr  | nSgnP | tSgnP | nSgnG | tSgnG | time\n");
    fclose(f);
 }
 else
 {
  printf("Syntax: chull <Dateiname>\n");
  exit(1);
 }
}