/* geo-ext3.c
 * 3D Geometric extentions to Postgres developed by TNO, The Netherlands.
 * In analogy with the 2D geometric data types developed by the
 * Postgres Research group, Berkeley.
 */

#include <stdio.h>
#include "utils/log.h"
#include "utils/geo-decls.h"
#include "geo-ext3.h"

#define LDELIM          '('
#define RDELIM          ')'
#define DELIM           ','
#define MDELIM          ':'
#define POINT3NARGS      3
#define LSEG3NARGS       6
#define BOX3NARGS        6

#define POLY3ALLOC(N) \
(POLY3STRUCT *) PALLOC((unsigned) (sizeof(POLY3STRUCT) + \
                            (((N)-1) > 0 ? ((N)-1) : 0) \
                            * sizeof(POINT3)))

POINT3 *
Point3_in(str)
   char     *str;
{
   extern double     atof();
   char     *coord[POINT3NARGS], *p;
   int      i;
   POINT3   *result;

   if (str == NULL) {
      elog(WARN, "Empty POINT3 string");
      return(NULL);
   }
   for (i = 0, p = str; *p && i < POINT3NARGS && *p != RDELIM; p++)
      if (*p == DELIM || (*p == LDELIM && !i))
         coord[i++] = p + 1;
   if (i < POINT3NARGS - 1) {
      elog(WARN, "Ill formatted POINT3 string: %s", str);
      return(NULL);
   }
   result = (POINT3 *) palloc(sizeof (POINT3));
   result->x = (float) atof(coord[0]);
   result->y = (float) atof(coord[1]);
   result->z = (float) atof(coord[2]);
printf("(%g,%g,%g)\n", result->x, result->y, result->z);
   return(result);
}     


char *
Point3_out(pt)
   POINT3   *pt;
{
   char    *result;
 
   if (pt == NULL) {
      elog(WARN, "NULL pointer for POINT3");
      return(NULL);
   }
   result = (char *) palloc(40);
   (void) sprintf(result, "(%g,%g,%g)", pt->x, pt->y, pt->z);
printf("(%g,%g,%g)", pt->x, pt->y, pt->z);
   return(result);
}


LSEG3  *
Lseg3_in(str)
   char    *str;
{
   char    *coord[LSEG3NARGS], *p;
   int     i;
   LSEG3   *result;
   extern  double   atof();


   if (str == NULL) {
      elog(WARN, "NULL string LSEG3");
      return(NULL);
   }
   for (i = 0, p = str; *p && i < LSEG3NARGS && *p != RDELIM; p++)
      if (*p == DELIM || (*p == LDELIM && !i))
         coord[i++] = p + 1;
   if (i < LSEG3NARGS - 1) {
      elog(WARN, "Ill formatted LSEG3: %s", str);
      return(NULL);
   }
   result = PALLOCTYPE(LSEG3);
   result->p[0].x = (float) atof(coord[0]);
   result->p[0].y = (float) atof(coord[1]);
   result->p[0].z = (float) atof(coord[2]);
   result->p[1].x = (float) atof(coord[3]);
   result->p[1].y = (float) atof(coord[4]);
   result->p[1].z = (float) atof(coord[5]);

   return(result);
}
 
char *
Lseg3_out(ls)
   LSEG3    *ls;
{
   char     *result; 
 
   if (ls == NULL) {
      elog(WARN, "NULL pointer for LSEG3");
      return(NULL);
   }
   result = (char *) palloc(80);
   (void) sprintf(result, "(%g,%g,%g,%g,%g,%g)",
      ls->p[0].x, ls->p[0].y, ls->p[0].z, ls->p[1].x, ls->p[1].y, ls->p[1].z);
   return(result);
}


BOX3 *
Box3_in(str)
   char    *str;
{
   char    *coord[BOX3NARGS], *p;
   int     i;
   BOX3   *result;
   extern  double   atof();
   float   tmp;

   if (str == NULL) {
      elog(WARN, "Empty string BOX3");
      return(NULL);
   }
   for (i = 0, p = str; *p && i < BOX3NARGS && *p != RDELIM; p++)
      if (*p == DELIM || (*p == LDELIM && !i))
         coord[i++] = p + 1;
   if (i < BOX3NARGS - 1) {
      elog(WARN, "Ill formatted BOX3: %s", str);
      return(NULL);
   }
   result = PALLOCTYPE(BOX3);
   result->xh = (float) atof(coord[0]);
   result->yh = (float) atof(coord[1]);
   result->zh = (float) atof(coord[2]);
   result->xl = (float) atof(coord[3]);
   result->yl = (float) atof(coord[4]);
   result->zl = (float) atof(coord[5]);

   if (result->xh < result->xl) {
      tmp = result->xh; result->xh = result->xl; result->xl = tmp;
   }
   if (result->yh < result->yl) {
      tmp = result->yh; result->yh = result->yl; result->yl = tmp;
   }
   if (result->zh < result->zl) {
      tmp = result->zh; result->zh = result->zl; result->zl = tmp;
   }

   return(result);
}
 
char *
Box3_out(box)
   BOX3    *box;
{
   char     *result; 
 
   if (box == NULL) {
      elog(WARN, "NULL pointer for BOX3");
      return(NULL);
   }
   result = (char *) palloc(80);
   (void) sprintf(result, "(%g,%g,%g,%g,%g,%g)",
      box->xl, box->yl, box->zl, box->xh, box->yh, box->zh);
   return(result);
}



/* Box3 functions*/
/*===============*/

POINT3 *
Center3Box(box)
   BOX3     *box;
{
   POINT3   *result;

   result = PALLOCTYPE(POINT3);
   result->x = (float)((box->xh + box->xl) / 2.0);
   result->y = (float)((box->yh + box->yl) / 2.0);
   result->z = (float)((box->zh + box->zl) / 2.0);

   return(result);
}

double *
Contents3Box(box)
   BOX3     *box;
{
   double   *result;

   result = PALLOCTYPE(double);

   *result = (double) 
           ((box->xh - box->xl) * (box->yh - box->yl) * (box->zh - box->zl));
   return(result);
}

BOX3 *
Intersect3BoxBox(box1, box2)
   BOX3     *box1, *box2;
{
   BOX3     *result;
   long     Overlap3BoxBox();

   if (! Overlap3BoxBox(box1,box2)) {
/*      elog(WARN, "No intersection BOX3 BOX3"); */
      return(NULL);
   }
   result = PALLOCTYPE(BOX3);
   result->xh = MIN(box1->xh, box2->xh);
   result->xl = MAX(box1->xl, box2->xl);
   result->yh = MIN(box1->yh, box2->yh);
   result->yl = MAX(box1->yl, box2->yl);
   result->zh = MIN(box1->zh, box2->zh);
   result->zl = MAX(box1->zl, box2->zl);

   return(result);
}

BOX3 *
Union3BoxBox(box1, box2)
   BOX3     *box1, *box2;
{
   BOX3     *result;

   result = PALLOCTYPE(BOX3);
   result->xh = MAX(box1->xh, box2->xh);
   result->xl = MIN(box1->xl, box2->xl);
   result->yh = MAX(box1->yh, box2->yh);
   result->yl = MIN(box1->yl, box2->yl);
   result->zh = MAX(box1->zh, box2->zh);
   result->zl = MIN(box1->zl, box2->zl);

   return(result);
}

long 
Equal3BoxBox(box1, box2)
   BOX3     *box1, *box2;
{
   return(FPeq((double)box1->xh, (double)box2->xh) &&
          FPeq((double)box1->xl, (double)box2->xl) &&
          FPeq((double)box1->yh, (double)box2->yh) &&
          FPeq((double)box1->yl, (double)box2->yl) &&
          FPeq((double)box1->zh, (double)box2->zh) &&
          FPeq((double)box1->zl, (double)box2->zl));
}

long
Overlap3BoxBox(box1, box2)
   BOX3     *box1, *box2;
{
   return(((box1->xh >= box2->xh && box1->xl <= box2->xh) ||
                (box2->xh >= box1->xh && box2->xl <= box1->xh)) &&
               ((box1->yh >= box2->yh && box1->yl <= box2->yh) ||
                (box2->yh >= box1->yh && box2->yl <= box1->yh)) &&
               ((box1->zh >= box2->zh && box1->zl <= box2->zh) ||
                (box2->zh >= box1->zh && box2->zl <= box1->zh)) );
}

long
Overleft3BoxBox(box1, box2)
   BOX3     *box1, *box2;
{
   return(box1->xh <= box2->xh);
}

long
Left3BoxBox(box1, box2)
   BOX3     *box1, *box2;
{
   return(box1->xh < box2->xh);
}

long
Overright3BoxBox(box1, box2)
   BOX3     *box1, *box2;
{
   return(box1->xh >= box2->xh);
}

long
Right3BoxBox(box1, box2)
   BOX3     *box1, *box2;
{
   return(box1->xh > box2->xh);
}

long
Contain3BoxBox(box1, box2)
   BOX3     *box1, *box2;
{
   return((box1->xh >= box2->xh && box1->xl <= box2->xl &&
           box1->yh >= box2->yh && box1->yl <= box2->yl &&
           box1->zh >= box2->zh && box1->zl <= box2->zl));
}

long
Contained3BoxBox(box1, box2)
   BOX3     *box1, *box2;
{
   return((box1->xh <= box2->xh && box1->xl >= box2->xl &&
           box1->yh <= box2->yh && box1->yl >= box2->yl &&
           box1->zh <= box2->zh && box1->zl >= box2->zl));
}

long
Above3BoxBox(box1, box2)
   BOX3     *box1, *box2;
{
   return(box1->yl >= box2->yh);
}

long
Below3BoxBox(box1, box2)
   BOX3     *box1, *box2;
{
   return(box1->yh <= box2->yl);
}

long
Higher3BoxBox(box1, box2)
   BOX3     *box1, *box2;
{
   return(box1->zl >= box2->zh);
}

long
Lower3BoxBox(box1, box2)
   BOX3     *box1, *box2;
{
   return(box1->zh <= box2->zl);
}

/* Point3 functions*/
/*=================*/

long
On3PntPgn(pt, box)
   POINT3   *pt;
   BOX3     *box;
{
   return( pt->x <= box->xh && pt->x >= box->xl &&
           pt->y <= box->yh && pt->y >= box->yl &&
           pt->z <= box->zh && pt->z >= box->zl );
}


long
Equal3PntPnt(pt1, pt2)
   POINT3   *pt1, *pt2;
{
   return( FPeq( (double)(pt1->x), (double)(pt2->x) ) &&
      FPeq( (double)(pt1->y), (double)(pt2->y) ) &&
      FPeq( (double)(pt1->z), (double)(pt2->z) ) );
}

 /* DISTANCE */
 /* ======== */
  
   
double
point3_dt(pt1, pt2)
   POINT3   *pt1, *pt2;
{
   return( sqrt( 
	  (double) (pt1->x - pt2->x)*(pt1->x - pt2->x) +
      (double) (pt1->y - pt2->y)*(pt1->y - pt2->y) +
	  (double) (pt1->z - pt2->z)*(pt1->z - pt2->z)));
}

double *
Distance3PntPnt(p1, p2)
   POINT3  *p1, *p2;
{
   double *result;
   double  point3_dt();

   result = PALLOCTYPE(double);
   *result = point3_dt(p1, p2);

   return( result );
}


/* POLY3 STRUCTURES */
/*==================*/


/*----------------------------------------------------------
 *  String to path3 / path3 to string conversion.
 *      External format:
 *              "(poly.npts: xcoord, ycoord, zcoord... )"
 *---------------------------------------------------------*/
       
POLY3STRUCT *
path3_in(str)
   char    *str;
{
   float   coord;
   extern double   atof();
   extern long     atol();
   long         field;
   char    *s;
   int     ct, i;
   POLY3STRUCT   *result;
 
   if (str == NULL) {
      elog(WARN, "Empty POLY3STRUCT string");
      return(NULL);
   }
   /* read the path3 header information */
   for (i = 0, s = str; *s && i < 1 && *s != RDELIM ; ++s)
      if (*s == DELIM || (*s == LDELIM && !i)) {
         field = atol(s + 1);
         ++i;
      }
   if (i < 0) {
      elog(WARN, "Error POLY3STRUCT string: %s", str);
      return(NULL);
   }
   result = POLY3ALLOC(field);
   result->npts =  field;

   /* read the path3 points */
 
   ct = result->npts * 3;  /* three coords for every point */
   for (i = 0; *s && i < ct && *s != RDELIM; ++s) {
      if (*s == DELIM || *s == MDELIM) {
         coord = (float) atof(s + 1);
         if ((i % 3) == 0)
            (result->p[i/3]).x = coord;
         else if ((i % 3) == 1)
            (result->p[i/3]).y = coord;
         else
            (result->p[i/3]).z = coord;
         ++i;
      }
   }
   if (i % 3 || i < --ct) {
      PFREE(result);
      elog(WARN, "Ill formatted POLY3STRUCT: %s", str);
      return(NULL);
   }
 
#ifdef DEBUG
   printf("path3_in:  %d: ", result->npts);
   for (i = 0; i < result->npts; ++i)
      printf("%g %g %g, ", result->p[i].x, result->p[i].y, result->p[i].z);
   putchar('\n');
#endif
 
   return(result);
}
 
 
char *
path3_out(pth)
   POLY3STRUCT    *pth;
{
   char            buf[50000], *result, *s;       
   int             i;

   static char     tmp[64];
 
   if (pth == NULL) {
      elog(WARN, "Empty string POLY3STRUCT");
      return(NULL);
   }
   (void) sprintf(buf,"%c%d\0", LDELIM, pth->npts);
   s = buf + strlen(buf);
   (void) sprintf(tmp, ":%g,%g,%g", pth->p[0].x, pth->p[0].y, pth->p[0].z);
   (void) strcpy(s, tmp);
   s += strlen(tmp);
   for (i = 1; i < pth->npts; ++i) {
      (void) sprintf(tmp, ",%g,%g,%g", pth->p[i].x, pth->p[i].y, pth->p[i].z);
      (void) strcpy(s, tmp);
      s += strlen(tmp);
   }
   *s++ = RDELIM;
   *s = '\0';
#ifdef DEBUG
   printf("path3_out: buf = %s\n", buf);
#endif
   result = (char *) palloc(strlen(buf) + 1);
   (void) strcpy(result, buf);

   return(result);
}

POLYLINE3 *
Polyline3_in(str)
   char *str;
{
   return((POLYLINE3 *) path3_in (str));
}

POLYGON3 *
Polygon3_in(str)
   char *str;
{
   return((POLYGON3 *) path3_in(str));
}

char *
Polyline3_out(pline)
   POLYLINE3  *pline;
{
   return(path3_out((POLY3STRUCT *) pline));
}

char *
Polygon3_out(pgon)
   POLYGON3  *pgon;
{
   return(path3_out((POLY3STRUCT *) pgon));
}


POINT3 *
point3_construct(x, y, z)
   double  x, y, z;
{
   POINT3   *result;

   result = PALLOCTYPE(POINT3);
   result->x = (float) (x);
   result->y = (float) (y);
   result->z = (float) (z);
   return(result);
}
 
 
long
On3PntPln(pt, pth)
   POINT3   *pt;
   POLYLINE3    *pth;
{
   int      i, n;
   double   point3_dt(), a, b;
 
   n = pth->npts - 1;
   a = point3_dt(pt, &pth->p[0]);
   for (i = 0; i < n; i++) {
      b = point3_dt(pt, &pth->p[i+1]);
      if (FPeq(a+b, point3_dt(&pth->p[i], &pth->p[i+1])))
         return(1);
      a = b;
   }
   return(0);
}

/* AREA */
/*======*/
 
double * 
area3(pgon)
   POLYGON3 *pgon;
{
/* to be implemented */
   return((double *)NULL);
} 


double *
chain3_length(pth)
   POLY3STRUCT *pth;
{  
   double   *result;
   int       i;

   result = PALLOCTYPE(double);
  *result = 0;
   for (i = 0; i < (pth->npts-1); ++i)
       *result += *Distance3PntPnt(&pth->p[i], &pth->p[i+1]);

   return(result);
}


POINT3 *
point3_copy(pt)
   POINT3   *pt;
{
   POINT3   *result;

   result = PALLOCTYPE(POINT3);
   result->x = pt->x;
   result->y = pt->y;
   result->z = pt->z;
   return(result);
}

double *
Length3Pln(pln)
   POLYLINE3  *pln;
{
   double     *chain3_length();

   return(chain3_length((POLY3STRUCT *) pln));
}              

long
Equal3PlnPln(pl1, pl2)
   POLYLINE3    *pl1, *pl2;
{
   long          n, result;
   long          Equal3PntPnt();
   int           i;
        
   if (pl1->npts == pl2->npts)
   {
      result = 1;
      n = pl1->npts;
      if (Equal3PntPnt(&pl1->p[0], &pl2->p[0]))
      {
         for (i = 0; i < n; i++)
            result = (result && Equal3PntPnt(&pl1->p[i], &pl2->p[i]));
      }
      else if (Equal3PntPnt(&pl1->p[0], &pl2->p[n-1]))
      {  
         for (i = 0; i < n; i++)
         result = (result && 
         Equal3PntPnt(&pl1->p[i], &pl2->p[n-1-i]));
      }
      else result = 0;
   }    
   else result = 0; 

   return (result);
}


POINT3 *
Point3Selector(pline, nr)
   POLYLINE3 *pline;
   long        nr;
{
   POINT3 *result;

   result = PALLOCTYPE(POINT3);

   if (0 < nr && nr <= pline->npts)
     nr--; 
   else if (nr < 0 && nr >= -pline->npts)
     nr = pline->npts + nr;
   else
     (void) printf("Not a legal number");
   result->x  = pline->p[nr].x;
   result->y  = pline->p[nr].y;
   result->z  = pline->p[nr].z;
   return(result);
}



double *
Area3Pgn(pgon)
   POLYGON3  *pgon;
{
/* to be implemented */
   return((double *)NULL);
}


double  *
Perimeter3Pgn(pgon)
   POLYGON3    *pgon;
{
   double     *Distance3PntPnt(), *result;
   double     *chain3_length();

   result = chain3_length((POLY3STRUCT *) pgon);
   *result += *Distance3PntPnt(&pgon->p[pgon->npts-1], &pgon->p[0]);
   return(result);
}


long
Equal3PgnPgn(pl1, pl2)
   POLYGON3    *pl1, *pl2;
{
   long         n, result;
   long         Equal3PntPnt();
   int          i, j;
 
   if (pl1->npts == pl2->npts)
   {
      result = 1;
      n = pl1->npts;
      i = 0;
      while (Equal3PntPnt(&pl1->p[0], &pl2->p[i]) && (i < n))
              i++;
      i--;
      if (Equal3PntPnt(&pl1->p[0], &pl2->p[i]))
      {
         if (Equal3PntPnt(&pl1->p[1], &pl2->p[i+1 % n]))
            for (j = 0; j < n; j++) 
               result = (result && Equal3PntPnt(&pl1->p[j], &pl2->p[i+j % n])); 
         else if (Equal3PntPnt(&pl1->p[1], &pl2->p[i-1 % n]))
            for (j = 0; j < n; j++) 
               result = (result && Equal3PntPnt(&pl1->p[j], &pl2->p[i-j % n]));
      }
   }
   else result = 0;

   return (result);
}