---

geoobj.c

Go to the documentation of this file.

/* geoobj.c -- A generic spatial object type for FMaps 
 *
 * Copyright (c) 2000 Eric G. Miller <egm2@jps.net>           -- for initial code
 *                    Franck Martin <franck@sopac.org>        
 *                    Gene Selkov Jr <selkovjr@mcs.anl.gov>   -- for GIST functions
 * 
 * This Program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This Program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU General Public
 * License along with this Program; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

/*
 * Modified by the FMaps Team and others 2000.  See the AUTHORS
 * file for a list of people on the FMaps Team.  See the ChangeLog
 * files for a list of changes.  These files are distributed with
 * FMaps at http://FMaps.sourceforge.net/.
 */

#include <stdlib.h>
#include <ctype.h>
#include <float.h>
#include <postgres.h>
#include <libpq-fe.h>

#include <access/gist.h>
#include <access/itup.h>
#include <access/rtree.h>

#include "geoobj.h"

/* To implement a new type modify the following procedures:
 *
 * getFMTypeCode 
 * getFMTypeStr
 * GeoObjParser
 * geoobj_out
 * Envelope
 */

/* Debugging */
#define TRACE(a) elog(NOTICE,a);

/* For buffer string */
#define MAX_FMTYPE_STRING 20

/* For the BIG/LITTLE ENDIAN */
#define BYTEORDER 1

/* For max double string length needs float.h */
#define P_MAXDIG DBL_DIG
#define P_MAXLEN (2*(DBL_DIG+7)+1)
static int  digits8 = DBL_DIG;

#define max(a,b)        ((a) >  (b) ? (a) : (b))
#define min(a,b)        ((a) <= (b) ? (a) : (b))
#define abs(a)          ((a) <  (0) ? (-a) : (a))

/* Forward Declarations for things not in geoobj.h */
int getFMTypeCode(char *begin, char **end, int4 *type);
int getFMTypeStr(GeoObj *obj, char *type);
Point *GeoObjPoints(char *begin, char** end, int4 *numpoints);
WKBGeometry *GeoObjParser(char *begin ,char **end, int4 type, int4 *size);
GeoObj *newGeoObj(char *instr);
int makeBox(char *input, char *output);

/* Input/Output for GeoObj */

GeoObj *geoobj_in (char *);
char *geoobj_out (GeoObj *);

/* Return the Bounding Box of a GeoObj as a GeoObj */
GeoObj *envelope (GeoObj *);

/* 
** GiST support methods
*/

bool             g_geoobj_consistent(GISTENTRY *entry, GeoObj *query, StrategyNumber strategy);
GISTENTRY *      g_geoobj_compress(GISTENTRY *entry);
GISTENTRY *      g_geoobj_decompress(GISTENTRY *entry);
float *          g_geoobj_penalty(GISTENTRY *origentry, GISTENTRY *newentry, float *result);
GIST_SPLITVEC *  g_geoobj_picksplit(bytea *entryvec, GIST_SPLITVEC *v);
bool             g_geoobj_leaf_consistent(GeoObj *key, GeoObj *query, StrategyNumber strategy);
bool             g_geoobj_internal_consistent(GeoObj *key, GeoObj *query, StrategyNumber strategy);
GeoObj *         g_geoobj_union(bytea *entryvec, int *sizep);
GeoObj *         g_geoobj_binary_union(GeoObj *r1, GeoObj *r2, int *sizep);
bool *           g_geoobj_same(GeoObj *b1, GeoObj *b2, bool *result);

/*
** R-tree suport functions
*/
bool         geoobj_same(GeoObj *a, GeoObj *b);
bool         geoobj_different(GeoObj *a, GeoObj *b);
bool         geoobj_contains(GeoObj *a, GeoObj *b);
bool         geoobj_contained (GeoObj *a, GeoObj *b);
bool         geoobj_overlap(GeoObj *a, GeoObj *b);
GeoObj *     geoobj_union(GeoObj *a, GeoObj *b);
GeoObj *     geoobj_inter(GeoObj *a, GeoObj *b);
float *      geoobj_size(GeoObj *a);
void         rt_geoobj_size(GeoObj *a, float *sz);

/*
** These make no sense for this type, but R-tree wants them
*/
bool         geoobj_over_left(GeoObj *a, GeoObj *b);
bool         geoobj_over_right(GeoObj *a, GeoObj *b);
bool         geoobj_left(GeoObj *a, GeoObj *b);
bool         geoobj_right(GeoObj *a, GeoObj *b);

/* 
** Auxiliary funxtions
*/
static       float distance_1D(float a1, float a2, float b1, float b2);
static       GeoObj *swap_corners (GeoObj *a);

int getFMTypeCode(char *begin, char **end, int4 *type)
{
  int i;
  char *cp, *p;
  char buffer[MAX_FMTYPE_STRING + 1];

TRACE("getFMTypeCode")

  p = begin;
  *type = 0; 
  while(isspace(*p))
    p++;
  cp = p;
  while(isalpha(*p))
    p++;
  if (p - cp > MAX_FMTYPE_STRING)
    return *type;
  for(i = 0; cp < p && i < MAX_FMTYPE_STRING; i++, cp++)
    buffer[i] = *cp;
  buffer[i] = '\0';
 
  if (strcasecmp(buffer,"POINT") == 0)
    *type = FM_POINT;
  else if (strcasecmp(buffer,"LINESTRING") == 0)
    *type = FM_LINESTRING;
  else if (strcasecmp(buffer,"POLYGON") == 0)
    *type = FM_POLYGON;
  else if (strcasecmp(buffer, "MULTIPOINT") == 0)
    *type = FM_MULTIPOINT;
  else if (strcasecmp(buffer,"MULTILINESTRING") == 0)
    *type = FM_MULTILINESTRING;
  else if (strcasecmp(buffer,"MULTIPOLYGON") == 0)
    *type = FM_MULTIPOLYGON;
  else if (strcasecmp(buffer,"BOX") == 0)
    *type = FM_BOX;   
    
  *end = cp;
  return *type;
}

/* getFMTypeStr: returns the string representation (uppercase ASCII)
 * for the GeoObj type, stored as an integer code in the structure.
 * On failure it returns zero.  Useful for a function to handle queries
 * that want the type string to filter unwanted objects.
 */

int getFMTypeStr(GeoObj *obj, char *type)
{

TRACE("getFMTypeStr")

  if(obj == NULL || type == NULL)
    return 0;

  switch(obj->wkbtype) {
  case FM_POINT: 
    strncpy(type, "POINT", MAX_FMTYPE_STRING);
    return FM_POINT;
  case FM_LINESTRING:
    strncpy(type, "LINESTRING", MAX_FMTYPE_STRING);
    return FM_LINESTRING;
  case FM_POLYGON:
    strncpy(type,"POLYGON", MAX_FMTYPE_STRING);
    return FM_POLYGON;
  case FM_MULTILINESTRING:
    strncpy(type,"MULTILINESTRING", MAX_FMTYPE_STRING);
    return FM_MULTILINESTRING;
  case FM_MULTIPOLYGON:
    strncpy(type,"MULTIPOLYGON", MAX_FMTYPE_STRING);
    return FM_MULTIPOLYGON;
  case FM_BOX:
    strncpy(type,"BOX", MAX_FMTYPE_STRING);
    return FM_BOX;
  default:
    strncpy(type,"", MAX_FMTYPE_STRING);
  }
  return 0;
}

/* GeoObjPoints: retreive an array of points */

Point *GeoObjPoints(char *begin, char** end, int4 *numpoints)
{
  char *p,*cp,*strend;
  char buffer[255];
  Point *points=NULL;
  double x,y,z;

TRACE("GeoObjPoints")


  strend = begin;
  while(*strend != '\0')
    strend++;
  
  p = begin;
  cp = p;  
  *numpoints=0;  
  

  while(isspace(*p))
    p++,cp++;
  
  if (*p != '(')
    elog(ERROR,"Syntax Error: Expecting (");
  p++,cp++;

  while(p < strend)
  {  
    while(isspace(*p))
      p++, cp++; 
    
    x = strtod(p,&cp);
    (p != cp) ? p = cp : elog(ERROR, "Numeric conversion error on first coordinate"); 
  
    while(isspace(*p))
      p++, cp++; 
         
    y = strtod(p,&cp);
    (p != cp) ? p = cp : elog(ERROR, "Numeric conversion error on second coordinate");  
  
    while(isspace(*p))
      p++, cp++; 
      
    z = strtod(p,&cp);
    (p != cp) ? p = cp : elog(ERROR, "Numeric conversion error on third coordinate");  
    
    *numpoints+=1;

    points=realloc(points,sizeof(Point)*(*numpoints));
    points[*numpoints-1].x = x;
    points[*numpoints-1].y = y;
    points[*numpoints-1].z = z;

sprintf(buffer,"points: %.g %.g %.g",points[*numpoints-1].x,points[*numpoints-1].y,points[*numpoints-1].z);
elog(NOTICE,buffer);

    while(isspace(*p))
      p++, cp++;
      
    if (*p != ',' && *p != ')')
      elog(ERROR,"Syntax Error: Expecting ) or ,");
    if (*p==')') 
    {
      p++,cp++;
      break;
    }
    p++,cp++;  


  }
  *end=cp;
  return points;
  
}
/* GeoObjParser: This function parses an input string and return
 * a WKBGeometry that will be stored in a geoobj structure
 */
 
WKBGeometry *GeoObjParser(char *begin ,char **end, int4 type, int4 *size)
{
  WKBGeometry *wkbgeo;
  Point *points;
  int4 i,numpoints;
  char *p,*cp;
  char buffer[255];

TRACE("GeoObjParser")
  
  *size=0;

  /* initialise the parameters */
  switch(type)
  {
    case FM_POINT:
      points = GeoObjPoints(begin,&p,&numpoints);  
      if (numpoints!=1) elog(ERROR,"Too many nodes found for type point");
      
      *size=sizeof(WKBPoint);      
      wkbgeo=palloc(*size);
      wkbgeo->point.byteOrder = BYTEORDER;  
      wkbgeo->point.wkbType = type;
      wkbgeo->point.point.x = points[0].x;
      wkbgeo->point.point.y = points[0].y;
      wkbgeo->point.point.z = points[0].z;
      free(points);

      break;
      
    case FM_LINESTRING:
      points = GeoObjPoints(begin,&p,&numpoints);
      if (numpoints<=1) elog(ERROR,"Not Enough nodes for type linestring");
      *size=sizeof(Point)*(numpoints)+sizeof(WKBLineString);
      wkbgeo=palloc(*size);
      wkbgeo->linestring.byteOrder = BYTEORDER;
      wkbgeo->linestring.wkbType = type;
      wkbgeo->linestring.numPoints = numpoints;
      for(i=0;i<numpoints;i++)
      {
        wkbgeo->linestring.points[i].x = points[i].x;
        wkbgeo->linestring.points[i].y = points[i].y;
        wkbgeo->linestring.points[i].z = points[i].z;         
      }   
      free(points);
      break;

    case FM_POLYGON:
      p = begin;
      cp = p;
      while(isspace(*p))
       p++,cp++;
  
      if (*p != '(')
        elog(ERROR,"Syntax Error: Expecting ( after Object type");
      p++,cp++;
      
      points = GeoObjPoints(p,&cp,&numpoints);
      (p != cp) ? p = cp : elog(ERROR, "Cannot extract the ring"); 

      while(isspace(*p))
       p++,cp++;
  
      if (*p != ',' && *p != ')')
        elog(ERROR,"Syntax Error: Expecting ) or , at the end of the ring");
      p++,cp++;

      if (numpoints<=2) elog(ERROR,"Not Enough nodes for type Polygon");
      *size=(sizeof(Point)*(numpoints)+sizeof(LinearRing))+sizeof(WKBPolygon);
      wkbgeo=palloc(*size);
      wkbgeo->polygon.byteOrder = BYTEORDER;
      wkbgeo->polygon.wkbType = type;
      wkbgeo->polygon.numRings = 1;
      wkbgeo->polygon.rings[0].numPoints = numpoints;
      
      for(i=0;i<numpoints;i++)
      {
        wkbgeo->polygon.rings[0].points[i].x = points[i].x;
        wkbgeo->polygon.rings[0].points[i].y = points[i].y;
        wkbgeo->polygon.rings[0].points[i].z = points[i].z;         
      }   
      free(points);
      break;

    case FM_BOX:
      points = GeoObjPoints(begin,&p,&numpoints);  
      if (numpoints!=2) elog(ERROR,"Two nodes only are required for type box");
      
      *size=sizeof(WKBBox);      
      wkbgeo=palloc(*size);
      wkbgeo->box.byteOrder = BYTEORDER; 
      wkbgeo->box.wkbType = type;
      wkbgeo->box.x[0] = points[0].x;
      wkbgeo->box.x[1] = points[0].y;
      wkbgeo->box.x[2] = points[0].z;
      wkbgeo->box.x[3] = points[1].x;
      wkbgeo->box.x[4] = points[1].y;
      wkbgeo->box.x[5] = points[1].z;
      free(points);

      break;
            
    default:
      elog(ERROR,"Unknown Data Type passed in GeoObjParser");
  } 
sprintf(buffer,"size: %d type: %d",*size,type);
elog(NOTICE,buffer); 
  *end=p;
  return wkbgeo; 
  
}


/* newGeoObj: This is the main() of functions for creating a new
 * GeoObj structure from an input string. It mostly calls other
 * functions to do the work. It needs to call a parser to do
 * the string copy to remove extra whitespace (TODO).
 */

GeoObj *newGeoObj(char *instr)
{
  int4 size, length, type;
  char *p;
  GeoObj *obj;
  WKBGeometry *wkbgeo;
  char buffer[255];

TRACE("newGeoObj")

  /* retreive the type */
  p = instr;
  if (!getFMTypeCode(instr,&p,&type))
    elog(ERROR,"Unknown Data Type!");
  /* parse the rest and compact into a WKBGeometry */
  instr = p;
  wkbgeo = GeoObjParser(instr,&p,type, &length);
  
  /* store the result into the geoobj */
sprintf(buffer,"length: %d type:%d ",length,type);
elog(NOTICE,buffer);
  size = offsetof(GeoObj, data[0]) + length + 1;
  obj = palloc(size);
  memset(obj,0,size);
  obj->size = size;
  obj->wkbtype = type;
  obj->byte_order = BYTEORDER; 
  /* move the wkbgeo structure into the geoobj */
  memcpy(obj->data,wkbgeo,length);
  
  /* free wkbgeo*/
  pfree(wkbgeo);
  
  return obj;
}   

/*------------------DEPRECIATED----------------------------*/    
/* makeBox: Called with a pointer to an input string
 * pointing on or just before the first '(' and a 
 * preallocated output string. It returns zero if
 * either of the pointers are NULL, otherwise it
 * returns the dimensions unless it raises an error.
 */

int makeBox(char *input, char *output)
{
  float8 xMin,xMax,yMin,yMax,zMin,zMax,tmp;
  int dim;
  char *p, *end, *tail;

TRACE("makeBox")

  if (input == NULL || output == NULL)
    return 0;

  p = end = tail = input;
  xMin = yMin = zMin = DBL_MAX;
  xMax = yMax = zMax = DBL_MIN;
  dim = 0;
  
  while(*end != '\0')
    end++;
  
  while(p < end) {
    while(isspace(*p))
      p++, tail++;
    if(*p != '(')
      elog(ERROR, "Unknown character where '(' should have been!");
    p++,tail++;
    tmp = strtod(p,&tail);
    (p != tail) ? p = tail : elog(ERROR, "Numeric conversion error!");
    (tmp < xMin) ? xMin = tmp : 1 ;
    (tmp > xMax) ? xMax = tmp : 1 ;
    while(isspace(*p))
      p++, tail++;
    if(*p != ',')
      elog(ERROR, "Unknown character where ',' should have been!");
    p++, tail++;
    tmp = strtod(p,&tail);
    (p != tail) ? p = tail : elog(ERROR, "Numeric conversion error!");
    (tmp < yMin) ? yMin = tmp : 1 ;
    (tmp > yMax) ? yMax = tmp : 1 ;
    while(isspace(*p))
      p++, tail++;
    if (dim < 2) { /* checks will catch an error later */
      if (*p == ')') {
        dim = 2;
        zMin = zMax = 0.0;
      }
      else {
        dim = 3;
      } 
    }
    if (dim == 3) {
      if(*p != ',')
        elog(ERROR, "Unknown character where ',' should have been!");
      p++, tail++;
      tmp = strtod(p, &tail);
      (p != tail) ? p = tail : elog(ERROR, "Numeric conversion error!");
      (tmp < zMin) ? zMin = tmp : 1 ;
      (tmp > zMax) ? zMax = tmp : 1 ;
      while(isspace(*p))
        p++, tail++;
    }
    if(*p != ')')
      elog (ERROR, "Unknown character where ')' should have been!");
    p++, tail++;
    while(isspace(*p))
      p++, tail++;
  }
  sprintf(output,"(%.*g,%.*g,%.*g) (%.*g,%.*g,%.*g)",
          digits8, xMin, digits8, yMin, digits8, zMin, 
          digits8, xMax, digits8, yMax, digits8, zMax
          );

  return dim;
}

/* geoobj_in: The main entry point for new data from postgresql */

GeoObj *geoobj_in(char *instr)
{
TRACE("geoobj_in")
     return newGeoObj(instr);
}


/* geoobj_out: The main output function for the data that postgresql
 * calls directly. 
 * We have to rebuild the string.
 */

char *geoobj_out (GeoObj *obj)
{
  char *outstr,*begin,*coord,*middle;
  size_t length;
  int4 type;
  int4 i;
  int4 maxlength;
  WKBGeometry *wkbgeo;
  char buffer[255];
  char *data;

TRACE("geoobj_out")

  data = obj->data;
  wkbgeo = (WKBGeometry *) data;  

sprintf(buffer,"out: %d ",obj->wkbtype);
elog(NOTICE,buffer);

  switch (obj->wkbtype)
  {
    case FM_POINT:
      maxlength=20*3+MAX_FMTYPE_STRING+4;
      outstr=palloc(maxlength);
      snprintf(outstr,maxlength,"POINT (%.g %.g %.g)",wkbgeo->point.point.x,wkbgeo->point.point.y,wkbgeo->point.point.z);
      break;

    case FM_LINESTRING:
      maxlength=wkbgeo->linestring.numPoints*20*3+MAX_FMTYPE_STRING+4;
      outstr=palloc(maxlength);
      begin=palloc(MAX_FMTYPE_STRING+3);
      coord=palloc(20*3);
      middle=palloc(maxlength);
      
      snprintf(begin,MAX_FMTYPE_STRING,"LINESTRING (");
      snprintf(middle,maxlength,"");

      for (i=0;i<wkbgeo->linestring.numPoints;i++)
      {
        snprintf(coord,maxlength,"%.g %.g %.g",wkbgeo->linestring.points[i].x,wkbgeo->linestring.points[i].y,wkbgeo->linestring.points[i].z);
        middle=strncat(middle,coord,maxlength);
        if (i<wkbgeo->linestring.numPoints-1)
        {
          middle=strncat(middle,", ",maxlength);
        } 
      }
      snprintf(outstr,maxlength,"%s%s)",begin,middle);
      pfree(begin);
      pfree(coord);
      pfree(middle);
      break;
      
    case FM_POLYGON:
      maxlength=wkbgeo->polygon.numRings*wkbgeo->polygon.rings[0].numPoints*20*3+MAX_FMTYPE_STRING+4;
      outstr=palloc(maxlength);
      begin=palloc(MAX_FMTYPE_STRING+3);
      coord=palloc(20*3);
      middle=palloc(maxlength);
      
      snprintf(begin,MAX_FMTYPE_STRING,"POLYGON ((");
      snprintf(middle,maxlength,"");

      for (i=0;i<wkbgeo->polygon.rings[0].numPoints;i++)
      {
        snprintf(coord,maxlength,"%.g %.g %.g",wkbgeo->polygon.rings[0].points[i].x,wkbgeo->polygon.rings[0].points[i].y,wkbgeo->polygon.rings[0].points[i].z);
        middle=strncat(middle,coord,maxlength);
        if (i<wkbgeo->polygon.rings[0].numPoints-1)
        {
          middle=strncat(middle,", ",maxlength);
        } 
      }
      snprintf(outstr,maxlength,"%s%s))",begin,middle);
      pfree(begin);
      pfree(coord);
      pfree(middle);
      break;

    case FM_BOX:
      maxlength=20*6+MAX_FMTYPE_STRING+4;
      outstr=palloc(maxlength);
      snprintf(outstr,maxlength,"BOX (%.g %.g %.g, %.g %.g %.g)",wkbgeo->box.x[0],wkbgeo->box.x[1],wkbgeo->box.x[2],wkbgeo->box.x[3],wkbgeo->box.x[4],wkbgeo->box.x[5]);
      break;
            
    default: 
      elog(ERROR,"Unknown type in geoobj_out");
  }
elog(NOTICE,outstr);  
  return outstr;
}

/* ----------------------DEPRECIATED --------------------------*/
/* box_out: Makes a box for any object as a new GeoObj, but
 * handles the case where the type is already a box specially.
 */

GeoObj *bounds (GeoObj *obj)
{
  char    *outstr, *cp;
  int     dims, length;

TRACE("bounds")

  if (obj == NULL)
    elog (ERROR,"box_out was passed a NULL object!");
  
  if (obj->wkbtype == FM_BOX)
    return obj;
  else
    dims = 3;

  length = dims * P_MAXLEN * 2 + sizeof("BOX ");
  outstr = palloc(length + 1);
  memset(outstr,0,length + 1);
  strcpy(outstr,"BOX ");
  cp = outstr;
  while(*cp != '\0')
    cp++;
  if (!makeBox(strchr(obj->data,'('),cp))
    elog (ERROR, "makeBox failed!\n");

  return newGeoObj(outstr);
}

/* Envelope: Returns the Bounding Box of any geoobj as a BOX type
 *
 */

GeoObj *envelope (GeoObj *obj)
{
  GeoObj *result;
  char strbox[20*6+MAX_FMTYPE_STRING+4];
  double minx;
  double miny;
  double minz;
  double maxx;
  double maxy;
  double maxz;
  int4 i;
  char *data;
  WKBGeometry *wkbgeo;
  char buffer[255];

TRACE("envelope")
  
sprintf(buffer,"envelope: %d size:%d",obj->wkbtype,obj->size);  
elog(NOTICE,buffer);
  data = obj->data;
  wkbgeo = (WKBGeometry *) data;
  
  switch (obj->wkbtype)
  {
    case FM_POINT:
elog(NOTICE,"point");
      minx=maxx=wkbgeo->point.point.x;
      miny=maxy=wkbgeo->point.point.y;
      minz=maxz=wkbgeo->point.point.z;
      break;

    case FM_LINESTRING:
      minx=maxx=wkbgeo->linestring.points[0].x;
      miny=maxy=wkbgeo->linestring.points[0].y;
      minz=maxz=wkbgeo->linestring.points[0].z;
      
      for (i=0;i<wkbgeo->linestring.numPoints;i++)
      {
        if (wkbgeo->linestring.points[i].x<minx) minx=wkbgeo->linestring.points[i].x;
        if (wkbgeo->linestring.points[i].y<miny) miny=wkbgeo->linestring.points[i].y;
        if (wkbgeo->linestring.points[i].z<minz) minz=wkbgeo->linestring.points[i].z;
        if (wkbgeo->linestring.points[i].x>maxx) maxx=wkbgeo->linestring.points[i].x;
        if (wkbgeo->linestring.points[i].y>maxy) maxy=wkbgeo->linestring.points[i].y;
        if (wkbgeo->linestring.points[i].z>maxz) maxz=wkbgeo->linestring.points[i].z;
      }      
      break;
      
    case FM_POLYGON:
      minx=maxx=wkbgeo->polygon.rings[0].points[0].x;
      miny=maxy=wkbgeo->polygon.rings[0].points[0].y;
      minz=maxz=wkbgeo->polygon.rings[0].points[0].z;
      
      for (i=0;i<wkbgeo->polygon.rings[0].numPoints;i++)
      {
        if (wkbgeo->polygon.rings[0].points[i].x<minx) minx=wkbgeo->polygon.rings[0].points[i].x;
        if (wkbgeo->polygon.rings[0].points[i].y<miny) miny=wkbgeo->polygon.rings[0].points[i].y;
        if (wkbgeo->polygon.rings[0].points[i].z<minz) minz=wkbgeo->polygon.rings[0].points[i].z;
        if (wkbgeo->polygon.rings[0].points[i].x>maxx) maxx=wkbgeo->polygon.rings[0].points[i].x;
        if (wkbgeo->polygon.rings[0].points[i].y>maxy) maxy=wkbgeo->polygon.rings[0].points[i].y;
        if (wkbgeo->polygon.rings[0].points[i].z>maxz) maxz=wkbgeo->polygon.rings[0].points[i].z;
      }
      break;  
      
    case FM_BOX:    
      minx=wkbgeo->box.x[0];
      miny=wkbgeo->box.x[1];
      minz=wkbgeo->box.x[2];
      maxx=wkbgeo->box.x[3];
      maxy=wkbgeo->box.x[4];
      maxz=wkbgeo->box.x[5];
      break;
            
    default:
      elog(ERROR,"Unknown type in envelope");
  }
  sprintf(strbox,"BOX (%.g %.g %.g, %.g %.g %.g)",minx,miny,minz,maxx,maxy,maxz);
elog(NOTICE,strbox);
  result=newGeoObj(strbox);
  return(result);
}
/*****************************************************************************
 *                         GiST functions
 *****************************************************************************/
/*
** The GiST Consistent method for boxes
** Should return false if for all data items x below entry,
** the predicate x op query == FALSE, where op is the oper
** corresponding to strategy in the pg_amop table.
*/

bool g_geoobj_consistent(GISTENTRY *entry, GeoObj *query, StrategyNumber strategy)
{
TRACE("g_geoobj_consistent")

    /*
    ** if entry is not leaf, use g_geobj_internal_consistent,
    ** else use g_geoobj_leaf_consistent
    */
    if (GIST_LEAF(entry))
      return(g_geoobj_leaf_consistent((GeoObj *)(entry->pred), query, strategy));
    else
      return(g_geoobj_internal_consistent((GeoObj *)(entry->pred), query, strategy));
}

/*
** GiST Compress and Decompress methods for geoobj
** Compress object into its bounding box
** decompres the bounding box into itself
*/

GISTENTRY *g_geoobj_compress(GISTENTRY *entry)
{
  char buffer[255];
  GeoObj *obj;
TRACE("g_geoobj_compress")
obj=(GeoObj *)entry->pred;
sprintf(buffer,"type: %d",obj->wkbtype);
elog(NOTICE,buffer);
elog(NOTICE,geoobj_out(obj));
    entry->pred=(char *)envelope((GeoObj *)entry->pred);

    return(entry);
}

GISTENTRY *g_geoobj_decompress(GISTENTRY *entry)
{
  char buffer[255];
  GeoObj *obj;
TRACE("g_geoobj_decompress")
obj=(GeoObj *)entry->pred;
sprintf(buffer,"type: %d",obj->wkbtype);
elog(NOTICE,buffer);
elog(NOTICE,geoobj_out(obj));
    return(entry);
}

/* 
** SUPPORT ROUTINES
*/
bool g_geoobj_leaf_consistent(GeoObj *key, GeoObj *query, StrategyNumber strategy)
{
  bool retval;

TRACE("g_geoobj_leaf_consistent")
  /*
  fprintf(stderr, "leaf_consistent, %d\n", strategy);
  */
  switch(strategy) {
  case RTLeftStrategyNumber:
    retval = (bool)geoobj_left(key, query);
    break;
  case RTOverLeftStrategyNumber:
    retval = (bool)geoobj_over_left(key,query);
    break;
  case RTOverlapStrategyNumber:
    retval = (bool)geoobj_overlap(key, query);
    break;
  case RTOverRightStrategyNumber:
    retval = (bool)geoobj_over_right(key, query);
    break;
  case RTRightStrategyNumber:
    retval = (bool)geoobj_right(key, query);
    break;
  case RTSameStrategyNumber:
    retval = (bool)geoobj_same(key, query);
    break;
  case RTContainsStrategyNumber:
    retval = (bool)geoobj_contains(key, query);
    break;
  case RTContainedByStrategyNumber:
    retval = (bool)geoobj_contained(key,query);
    break;
  default:
    retval = FALSE;
  }
  return(retval);
}

bool g_geoobj_internal_consistent(GeoObj *key, GeoObj *query, StrategyNumber strategy)
{
  bool retval;

TRACE("g_geoobj_internal_consistent")
  
  /*
  fprintf(stderr, "internal_consistent, %d\n", strategy);
  */
  switch(strategy) {
  case RTLeftStrategyNumber:
  case RTOverLeftStrategyNumber:
    retval = (bool)geoobj_over_left(key,query);
    break;
  case RTOverlapStrategyNumber:
    retval = (bool)geoobj_overlap(key, query);
    break;
  case RTOverRightStrategyNumber:
  case RTRightStrategyNumber:
    retval = (bool)geoobj_right(key, query);
    break;
  case RTSameStrategyNumber:
  case RTContainsStrategyNumber:
    retval = (bool)geoobj_contains(key, query);
    break;
  case RTContainedByStrategyNumber:
    retval = (bool)geoobj_overlap(key, query);
    break;
  default:
    retval = FALSE;
    }
  return(retval);
}
/*
** End SUPPORT ROUTINES
*/


/*
** The GiST Union method for boxes
** returns the minimal bounding box that encloses all the entries in entryvec
*/
GeoObj *g_geoobj_union(bytea *entryvec, int *sizep)
{
    int numranges, i;
    GeoObj *out = (GeoObj *)NULL;
    GeoObj *tmp;
TRACE("g_geoobj_union")
    /*
    fprintf(stderr, "union\n");
    */
    numranges = (VARSIZE(entryvec) - VARHDRSZ)/sizeof(GISTENTRY); 
    tmp = (GeoObj *)(((GISTENTRY *)(VARDATA(entryvec)))[0]).pred;
    /*
    *sizep = sizeof(NDBOX); -- NDBOX has variable size
    */
    *sizep = tmp->size;

    for (i = 1; i < numranges; i++) {
      out = g_geoobj_binary_union(tmp, (GeoObj *)
                                 (((GISTENTRY *)(VARDATA(entryvec)))[i]).pred,
                                 sizep);
      /*
        fprintf(stderr, "\t%s ^ %s -> %s\n", cube_out(tmp), cube_out((NDBOX *)(((GISTENTRY *)(VARDATA(entryvec)))[i]).pred), cube_out(out));
      */
      if (i > 1) pfree(tmp);
      tmp = out;
    }

    return(out);
}

/*
** The GiST Penalty method for boxes
** As in the R-tree paper, we use change in area as our penalty metric
*/
float *g_geoobj_penalty(GISTENTRY *origentry, GISTENTRY *newentry, float *result)
{
    Datum ud;
    float tmp1, tmp2;

TRACE("g_geoobj_penalty")
    
    ud = (Datum)geoobj_union((GeoObj *)(origentry->pred), (GeoObj *)(newentry->pred));
    rt_geoobj_size((GeoObj *)ud, &tmp1);
    rt_geoobj_size((GeoObj *)(origentry->pred), &tmp2);
    *result = tmp1 - tmp2;
    pfree((char *)ud);
    /*
    fprintf(stderr, "penalty\n");
    fprintf(stderr, "\t%g\n", *result);
    */
    return(result);
}



/*
** The GiST PickSplit method for boxes
** We use Guttman's poly time split algorithm 
*/
GIST_SPLITVEC *g_geoobj_picksplit(bytea *entryvec, GIST_SPLITVEC *v)
{
    OffsetNumber i, j;
    GeoObj *datum_alpha, *datum_beta;
    GeoObj *datum_l, *datum_r;
    GeoObj *union_d, *union_dl, *union_dr;
    GeoObj *inter_d;
    bool firsttime;
    float size_alpha, size_beta, size_union, size_inter;
    float size_waste, waste;
    float size_l, size_r;
    int nbytes;
    OffsetNumber seed_1 = 0, seed_2 = 0;
    OffsetNumber *left, *right;
    OffsetNumber maxoff;

TRACE("g_geoobj_picksplit")

    /*
    fprintf(stderr, "picksplit\n");
    */
    maxoff = ((VARSIZE(entryvec) - VARHDRSZ)/sizeof(GISTENTRY)) - 2;
    nbytes =  (maxoff + 2) * sizeof(OffsetNumber);
    v->spl_left = (OffsetNumber *) palloc(nbytes);
    v->spl_right = (OffsetNumber *) palloc(nbytes);
    
    firsttime = true;
    waste = 0.0;
    
    for (i = FirstOffsetNumber; i < maxoff; i = OffsetNumberNext(i)) 
    {
            datum_alpha = (GeoObj *)(((GISTENTRY *)(VARDATA(entryvec)))[i].pred);
            for (j = OffsetNumberNext(i); j <= maxoff; j = OffsetNumberNext(j)) 
            {
              datum_beta = (GeoObj *)(((GISTENTRY *)(VARDATA(entryvec)))[j].pred);
            
              /* compute the wasted space by unioning these guys */
            /* size_waste = size_union - size_inter; */
            union_d = (GeoObj *)geoobj_union(datum_alpha, datum_beta);
            rt_geoobj_size(union_d, &size_union);
            inter_d = (GeoObj *)geoobj_inter(datum_alpha, datum_beta);
            rt_geoobj_size(inter_d, &size_inter);
            size_waste = size_union - size_inter;
            
            pfree(union_d);
            
            if (inter_d != (GeoObj *) NULL)
                pfree(inter_d);
            
              /*
             *  are these a more promising split than what we've
             *  already seen?
             */
            
              if (size_waste > waste || firsttime) 
              {
                waste = size_waste;
          seed_1 = i;
                      seed_2 = j;
                      firsttime = false;
              }
            }
    }
    
    left = v->spl_left;
    v->spl_nleft = 0;
    right = v->spl_right;
    v->spl_nright = 0;
    
    datum_alpha = (GeoObj *)(((GISTENTRY *)(VARDATA(entryvec)))[seed_1].pred);
    datum_l = (GeoObj *)geoobj_union(datum_alpha, datum_alpha);
    rt_geoobj_size((GeoObj *)datum_l, &size_l);
    datum_beta = (GeoObj *)(((GISTENTRY *)(VARDATA(entryvec)))[seed_2].pred);;
    datum_r = (GeoObj *)geoobj_union(datum_beta, datum_beta);
    rt_geoobj_size((GeoObj *)datum_r, &size_r);
    
    /*
     *  Now split up the regions between the two seeds.  An important
     *  property of this split algorithm is that the split vector v
     *  has the indices of items to be split in order in its left and
     *  right vectors.  We exploit this property by doing a merge in
     *  the code that actually splits the page.
     *
     *  For efficiency, we also place the new index tuple in this loop.
     *  This is handled at the very end, when we have placed all the
     *  existing tuples and i == maxoff + 1.
     */
    
    maxoff = OffsetNumberNext(maxoff);
    for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i)) 
    {
        
            /*
         *  If we've already decided where to place this item, just
         *  put it on the right list.  Otherwise, we need to figure
         *  out which page needs the least enlargement in order to
         *  store the item.
         */
        
            if (i == seed_1) 
            {
              *left++ = i;
              v->spl_nleft++;
              continue;
            } else if (i == seed_2) 
            {
              *right++ = i;
              v->spl_nright++;
              continue;
            }
        
            /* okay, which page needs least enlargement? */ 
        datum_alpha = (GeoObj *)(((GISTENTRY *)(VARDATA(entryvec)))[i].pred);
        union_dl = (GeoObj *)geoobj_union(datum_l, datum_alpha);
        union_dr = (GeoObj *)geoobj_union(datum_r, datum_alpha);
        rt_geoobj_size((GeoObj *)union_dl, &size_alpha);
        rt_geoobj_size((GeoObj *)union_dr, &size_beta);
        
        /* pick which page to add it to */
        if (size_alpha - size_l < size_beta - size_r) 
        {
              pfree(datum_l);
            pfree(union_dr);
            datum_l = union_dl;
            size_l = size_alpha;
            *left++ = i;
            v->spl_nleft++;
        } else 
        {
              pfree(datum_r);
            pfree(union_dl);
            datum_r = union_dr;
            size_r = size_alpha;
            *right++ = i;
            v->spl_nright++;
            }
    }
    *left = *right = FirstOffsetNumber; /* sentinel value, see dosplit() */
    
    v->spl_ldatum = (char *)datum_l;
    v->spl_rdatum = (char *)datum_r;

    return v;
}

/*
** Equality method
*/
bool *g_geoobj_same(GeoObj *b1, GeoObj *b2, bool *result)
{
  TRACE("g_geoobj_same")
  if (geoobj_same(b1, b2))
    *result = TRUE;
  else *result = FALSE;
  /*
  fprintf(stderr, "same: %s\n", (*result ? "TRUE" : "FALSE" ));
  */
  return(result);
}

GeoObj *g_geoobj_binary_union(GeoObj *r1, GeoObj *r2, int *sizep)
{
    GeoObj *retval;

TRACE("g_geoobj_binary_union")
    retval = geoobj_union(r1, r2);
    *sizep = retval->size;

    return (retval);
}


/* cube_union */
GeoObj *geoobj_union(GeoObj *box_a, GeoObj *box_b)
{
  int i;
  GeoObj *result;
  GeoObj *a = swap_corners(box_a);
  GeoObj *b = swap_corners(box_b);
  char *data;
  WKBGeometry *wkbgeoresult;
  WKBGeometry *wkbgeoa;
  WKBGeometry *wkbgeob;

TRACE("geoobj_union")  
  data=a->data;
  wkbgeoa=(WKBGeometry *)data;

  data=b->data;
  wkbgeob=(WKBGeometry *)data;
  
  result = palloc(a->size);
  result->size = a->size;
  result->byte_order = BYTEORDER;
  result->wkbtype = FM_BOX;
  data=result->data;
  wkbgeoresult=(WKBGeometry *)data;
  wkbgeoresult->box.byteOrder = BYTEORDER;
  wkbgeoresult->box.wkbType = FM_BOX;

  
  /* use the potentially smaller of the two boxes (b) to fill in 
     the result, padding absent dimensions with zeroes*/
  for ( i = 0; i < 3; i++ ) {
    wkbgeoresult->box.x[i] = wkbgeob->box.x[i];
    wkbgeoresult->box.x[i + 3] = wkbgeob->box.x[i + 3];
  }

    
  /* compute the union */
  for ( i = 0; i < 3; i++ ) {
    wkbgeoresult->box.x[i] = min(wkbgeoa->box.x[i], wkbgeoresult->box.x[i]);
  }
  for ( i = 3; i < 6; i++ ) {
    wkbgeoresult->box.x[i] = max(wkbgeoa->box.x[i], wkbgeoresult->box.x[i]);
  }

  pfree(a);
  pfree(b);

  return(result);
}

/* cube_inter */
GeoObj *geoobj_inter(GeoObj *box_a, GeoObj *box_b)
{
  int i;
  GeoObj * result;
  GeoObj *a = swap_corners(box_a);
  GeoObj *b = swap_corners(box_b);
  char *data;
  WKBGeometry *wkbgeoresult;
  WKBGeometry *wkbgeoa;
  WKBGeometry *wkbgeob;

TRACE("geoobj_inter")  
  data=a->data;
  wkbgeoa=(WKBGeometry *)data;

  data=b->data;
  wkbgeob=(WKBGeometry *)data;
    
  result = palloc(a->size);
  result->size = a->size;
  result->byte_order = BYTEORDER;
  result->wkbtype = FM_BOX;
  data=result->data;
  wkbgeoresult=(WKBGeometry *)data;
  wkbgeoresult->box.byteOrder = BYTEORDER;
  wkbgeoresult->box.wkbType = FM_BOX;

  /* use the potentially  smaller of the two boxes (b) to fill in 
     the result, padding absent dimensions with zeroes*/
  for ( i = 0; i < 3; i++ ) {
    wkbgeoresult->box.x[i] = wkbgeob->box.x[i];
    wkbgeoresult->box.x[i + 3] = wkbgeob->box.x[i + 3];
  }
    
  /* compute the intersection */
  for ( i = 0; i < 3; i++ ) {
    wkbgeoresult->box.x[i] = max(wkbgeoa->box.x[i], wkbgeoresult->box.x[i]);
  }
  for ( i = 3; i < 6; i++ ) {
    wkbgeoresult->box.x[i] = min(wkbgeoa->box.x[i], wkbgeoresult->box.x[i]);
  }
  
  pfree(a);
  pfree(b);

  /* Is it OK to return a non-null intersection for non-overlapping boxes? */
  return(result);
}

/* cube_size */
float *geoobj_size(GeoObj *a)
{
  int i,j;
  float *result;
  char *data;
  WKBGeometry *wkbgeoa;

TRACE("geoobj_size")  
  data=a->data;
  wkbgeoa=(WKBGeometry *)data;
  
  result = (float *) palloc(sizeof(float));
  
  *result = 1.0;
  for ( i = 0, j = 3; i < 3; i++,j++ ) {
    *result=(*result)*abs((wkbgeoa->box.x[j] - wkbgeoa->box.x[i]));
  }
  
  return(result);
}

void rt_geoobj_size(GeoObj *a, float *size)
{
  int i,j;
  char *data;
  WKBGeometry *wkbgeoa;

TRACE("rt_geoobj_size")  
  data=a->data;
  wkbgeoa=(WKBGeometry *)data;
  
  if (a == (GeoObj *) NULL)
    *size = 0.0;
  else {
    *size = 1.0;
    for ( i = 0, j = 3; i < 3; i++,j++ ) {
      *size=(*size)*abs((wkbgeoa->box.x[j] - wkbgeoa->box.x[i]));
    }
  }
  return;
}

/* The following four methods compare the projections of the boxes
   onto the 0-th coordinate axis. These methods are useless for dimensions
   larger than 2, but it seems that R-tree requires all its strategies
   map to real functions that return something */

/*  is the right edge of (a) located to the left of 
    the right edge of (b)? */
bool geoobj_over_left(GeoObj *box_a, GeoObj *box_b)
{
  GeoObj *a;
  GeoObj *b;
  char *data;
  WKBGeometry *wkbgeoa;
  WKBGeometry *wkbgeob;

TRACE("geoobj_over_left")  
  if ( (box_a==NULL) || (box_b==NULL) )
    return(FALSE);

  a = swap_corners(box_a);
  b = swap_corners(box_b);

  data=a->data;
  wkbgeoa=(WKBGeometry *)data;

  data=b->data;
  wkbgeob=(WKBGeometry *)data;
  
  return( wkbgeoa->box.x[3 - 1] <= wkbgeob->box.x[3 - 1] && !geoobj_left(a, b) && !geoobj_right(a, b) );
}

/*  is the left edge of (a) located to the right of 
    the left edge of (b)? */
bool geoobj_over_right(GeoObj *box_a, GeoObj *box_b)
{
  GeoObj *a;
  GeoObj *b;
  char *data;
  WKBGeometry *wkbgeoa;
  WKBGeometry *wkbgeob;

TRACE("geoobj_over_right")  
  if ( (box_a==NULL) || (box_b==NULL) )
    return(FALSE);

  a = swap_corners(box_a);
  b = swap_corners(box_b);

  data=a->data;
  wkbgeoa=(WKBGeometry *)data;

  data=b->data;
  wkbgeob=(WKBGeometry *)data; 

  return( wkbgeoa->box.x[3 - 1] >= wkbgeob->box.x[3 - 1] && !geoobj_left(a, b) && !geoobj_right(a, b) );
}


/* return 'true' if the projection of 'a' is
   entirely on the left of the projection of 'b' */
bool geoobj_left(GeoObj *box_a, GeoObj *box_b)
{
  GeoObj *a;
  GeoObj *b;
  char *data;
  WKBGeometry *wkbgeoa;
  WKBGeometry *wkbgeob;

TRACE("geoobj_left")    
  if ( (box_a==NULL) || (box_b==NULL) )
    return(FALSE);

  a = swap_corners(box_a);
  b = swap_corners(box_b);
  
  data=a->data;
  wkbgeoa=(WKBGeometry *)data;

  data=b->data;
  wkbgeob=(WKBGeometry *)data;
  
  return( wkbgeoa->box.x[3 - 1] < wkbgeob->box.x[0]);
}

/* return 'true' if the projection of 'a' is
   entirely on the right  of the projection of 'b' */
bool geoobj_right(GeoObj *box_a, GeoObj *box_b)
{
  GeoObj *a;
  GeoObj *b;
  char *data;
  WKBGeometry *wkbgeoa;
  WKBGeometry *wkbgeob;

TRACE("geoobj_right")    
  if ( (box_a==NULL) || (box_b==NULL) )
    return(FALSE);

  a = swap_corners(box_a);
  b = swap_corners(box_b);

  data=a->data;
  wkbgeoa=(WKBGeometry *)data;

  data=b->data;
  wkbgeob=(WKBGeometry *)data;

  return( wkbgeoa->box.x[0] > wkbgeob->box.x[3 - 1]);
}

/* make up a metric in which one box will be 'lower' than the other
   -- this can be useful for srting and to determine uniqueness */
bool geoobj_lt(GeoObj *box_a, GeoObj *box_b)
{
  int i;
  GeoObj *a;
  GeoObj *b;
  char *data;
  WKBGeometry *wkbgeoa;
  WKBGeometry *wkbgeob;

TRACE("geoobj_lt")    
  if ( (box_a==NULL) || (box_b==NULL) )
    return(FALSE);

  a = swap_corners(box_a);
  b = swap_corners(box_b);

  data=a->data;
  wkbgeoa=(WKBGeometry *)data;

  data=b->data;
  wkbgeob=(WKBGeometry *)data;

  /* if all common dimensions are equal, the cube with more dimensions wins */
  if ( geoobj_same(a, b) ) {
    return(FALSE);
  }

  /* compare the common dimensions */
  for ( i = 0; i < 3; i++ ) {
    if ( wkbgeoa->box.x[i] > wkbgeob->box.x[i] )
      return(FALSE); 
    if ( wkbgeoa->box.x[i] < wkbgeob->box.x[i] )
      return(TRUE); 
  }
  for ( i = 0; i < 3; i++ ) {
    if ( wkbgeoa->box.x[i + 3] > wkbgeob->box.x[i + 3] )
      return(FALSE); 
    if ( wkbgeoa->box.x[i +3] < wkbgeob->box.x[i + 3] )
      return(TRUE); 
  }

  return(FALSE);
}


bool geoobj_gt(GeoObj *box_a, GeoObj *box_b)
{
  int i;
  GeoObj *a;
  GeoObj *b;
  char *data;
  WKBGeometry *wkbgeoa;
  WKBGeometry *wkbgeob;

TRACE("geoobj_gt")    
  if ( (box_a==NULL) || (box_b==NULL) )
    return(FALSE);

  a = swap_corners(box_a);
  b = swap_corners(box_b);

  data=a->data;
  wkbgeoa=(WKBGeometry *)data;

  data=b->data;
  wkbgeob=(WKBGeometry *)data;
  
  /* if all common dimensions are equal, the cube with more dimensions wins */
  if ( geoobj_same(a, b) ) 
  {
    return(FALSE);
  }

  /* compare the common dimensions */
  for ( i = 0; i < 3; i++ ) 
  {
    if ( wkbgeoa->box.x[i] < wkbgeob->box.x[i] )
      return(FALSE); 
    if ( wkbgeoa->box.x[i] > wkbgeob->box.x[i] )
      return(TRUE); 
  }
  for ( i = 0; i < 3; i++ ) 
  {
    if ( wkbgeoa->box.x[i + 3] < wkbgeob->box.x[i + 3] )
      return(FALSE); 
    if ( wkbgeoa->box.x[i + 3] > wkbgeob->box.x[i + 3] )
      return(TRUE); 
  }

  return(FALSE);
}


/* Equal */
bool geoobj_same(GeoObj *box_a, GeoObj *box_b)
{
  int i;
  GeoObj *a;
  GeoObj *b;
  char *data;
  WKBGeometry *wkbgeoa;
  WKBGeometry *wkbgeob;

TRACE("geoobj_same")    
  if ( (box_a==NULL) || (box_b==NULL) )
    return(FALSE);

  a = swap_corners(box_a);
  b = swap_corners(box_b);
  
  data=a->data;
  wkbgeoa=(WKBGeometry *)data;

  data=b->data;
  wkbgeob=(WKBGeometry *)data;

  for ( i = 0; i < 3; i++ ) {
    if ( wkbgeoa->box.x[i] != wkbgeob->box.x[i] )
      return(FALSE);
    if ( wkbgeoa->box.x[i + 3] != wkbgeob->box.x[i + 3] )
      return(FALSE);
  }

  pfree(a);
  pfree(b);

  return(TRUE);
}

/* Different */
bool geoobj_different(GeoObj *box_a, GeoObj *box_b)
{
TRACE("geoobj_different")
  return(!geoobj_same(box_a, box_b));
}


/* Contains */
/* Box(A) CONTAINS Box(B) IFF pt(A) < pt(B) */
bool geoobj_contains(GeoObj *box_a, GeoObj *box_b)
{
  int i;
  GeoObj *a;
  GeoObj *b;
  char *data;
  WKBGeometry *wkbgeoa;
  WKBGeometry *wkbgeob;

TRACE("geoobj_contains")  
  if ( (box_a==NULL) || (box_b==NULL) )
    return(FALSE);

  a = swap_corners(box_a);
  b = swap_corners(box_b);

  data=a->data;
  wkbgeoa=(WKBGeometry *)data;

  data=b->data;
  wkbgeob=(WKBGeometry *)data;
  
  /* Can't care less about the excess dimensions of (a), if any */
  for ( i = 0; i < 3; i++ ) {
    if ( wkbgeoa->box.x[i] > wkbgeob->box.x[i] )
      return(FALSE);
    if ( wkbgeoa->box.x[i + 3] < wkbgeob->box.x[i + 3] )
      return(FALSE);
  }

  pfree(a);
  pfree(b);

  return(TRUE);
}

/* Contained */
/* Box(A) Contained by Box(B) IFF Box(B) Contains Box(A) */
bool geoobj_contained(GeoObj *a, GeoObj *b)
{
TRACE("geoobj_contained")
  if (geoobj_contains(b,a) == TRUE)
    return(TRUE);
  else
    return(FALSE);
}

/* Overlap */
/* Box(A) Overlap Box(B) IFF (pt(a)LL < pt(B)UR) && (pt(b)LL < pt(a)UR) */
bool geoobj_overlap(GeoObj *box_a, GeoObj *box_b)
{
  int i;
  GeoObj *a;
  GeoObj *b;
  char *data;
  WKBGeometry *wkbgeoa;
  WKBGeometry *wkbgeob;

TRACE("geoobj_overlap") 
elog(NOTICE,geoobj_out(box_a));
elog(NOTICE,geoobj_out(box_b));
 
  if ( (box_a==NULL) || (box_b==NULL) )
    return(FALSE);

/*  a = swap_corners(box_a);
  b = swap_corners(box_b);
*/

  a=envelope(box_a);
  b=envelope(box_b);

  data=a->data;
  wkbgeoa=(WKBGeometry *)data;

  data=b->data;
  wkbgeob=(WKBGeometry *)data;
  
  /* compare within the dimensions of (b) */
  for ( i = 0; i < 3; i++ ) {
    if ( wkbgeoa->box.x[i] > wkbgeob->box.x[i + 3] )
      return(FALSE);
    if ( wkbgeoa->box.x[i + 3] < wkbgeob->box.x[i] )
      return(FALSE);
  }

  pfree(a);
  pfree(b);

  return(TRUE);
}


/* Distance */
/* The distance is computed as a per axis sum of the squared distances
   between 1D projections of the boxes onto Cartesian axes. Assuming zero 
   distance between overlapping projections, this metric coincides with the 
   "common sense" geometric distance */
float *geoobj_distance(GeoObj *a, GeoObj *b)
{
  int i;
  double d, distance;
  float *result;
  char *data;
  WKBGeometry *wkbgeoa;
  WKBGeometry *wkbgeob;

TRACE("geoobj_distance")
  data=a->data;
  wkbgeoa=(WKBGeometry *)data;

  data=b->data;
  wkbgeob=(WKBGeometry *)data;
  
  result = (float *) palloc(sizeof(float));
  
  distance = 0.0;
  /* compute within the dimensions of (b) */
  for ( i = 0; i < 3; i++ ) {
    d = distance_1D(wkbgeoa->box.x[i], wkbgeoa->box.x[i + 3], wkbgeob->box.x[i], wkbgeob->box.x[i + 3]);
    distance += d*d;
  }

  *result = (float)sqrt(distance);

  return(result);
}

static float distance_1D(float a1, float a2, float b1, float b2)
{
TRACE("distance_1D")
  /* interval (a) is entirely on the left of (b) */
  if( (a1 <= b1) && (a2 <= b1) && (a1 <= b2) && (a2 <= b2) ) {
    return ( min( b1, b2 ) - max( a1, a2 ) );
  }

  /* interval (a) is entirely on the right of (b) */
  if( (a1 > b1) && (a2 > b1) && (a1 > b2) && (a2 > b2) ) {
    return ( min( a1, a2 ) - max( b1, b2 ) );
  }
  
  /* the rest are all sorts of intersections */
  return(0.0);
}

/* normalize the box's co-ordinates by placing min(xLL,xUR) to LL
   and max(xLL,xUR) to UR 
*/
static GeoObj *swap_corners( GeoObj *a )
{
  int i, j;
  GeoObj * result;
  char *data;
  WKBGeometry *wkbgeoa;
  WKBGeometry *wkbgeoresult;

TRACE("swap_corners")
  data=a->data;
  wkbgeoa=(WKBGeometry *)data;
    
  result = palloc(a->size);
  result->size = a->size;
  result->byte_order = BYTEORDER;
  result->wkbtype = FM_BOX;
  data=a->data;
  wkbgeoresult=(WKBGeometry *)data;
  wkbgeoresult->box.byteOrder = BYTEORDER;
  wkbgeoresult->box.wkbType = FM_BOX;
  
  for ( i = 0, j = 3; i < 3; i++, j++ ) {
    wkbgeoresult->box.x[i] = min(wkbgeoa->box.x[i],wkbgeoa->box.x[j]);
    wkbgeoresult->box.x[j] = max(wkbgeoa->box.x[i],wkbgeoa->box.x[j]);
  }

  return(result);
}

---