Main Page   Alphabetical List   Compound List   File List   Compound Members   File Members  

geoobj.c

Go to the documentation of this file.
00001 /* geoobj.c -- A generic spatial object type for FMaps 
00002  *
00003  * Copyright (c) 2000 Eric G. Miller <egm2@jps.net>           -- for initial code
00004  *                    Franck Martin <franck@sopac.org>        
00005  *                    Gene Selkov Jr <selkovjr@mcs.anl.gov>   -- for GIST functions
00006  * 
00007  * This Program is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU General Public
00009  * License as published by the Free Software Foundation; either
00010  * version 2 of the License, or (at your option) any later version.
00011  *
00012  * This Program is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00015  * Library General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU General Public
00018  * License along with this Program; if not, write to the
00019  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
00020  * Boston, MA 02111-1307, USA.
00021  */
00022 
00023 /*
00024  * Modified by the FMaps Team and others 2000.  See the AUTHORS
00025  * file for a list of people on the FMaps Team.  See the ChangeLog
00026  * files for a list of changes.  These files are distributed with
00027  * FMaps at http://FMaps.sourceforge.net/.
00028  */
00029 
00030 #include <stdlib.h>
00031 #include <ctype.h>
00032 #include <float.h>
00033 #include <postgres.h>
00034 #include <libpq-fe.h>
00035 
00036 #include <access/gist.h>
00037 #include <access/itup.h>
00038 #include <access/rtree.h>
00039 
00040 #include "geoobj.h"
00041 
00042 /* To implement a new type modify the following procedures:
00043  *
00044  * getFMTypeCode 
00045  * getFMTypeStr
00046  * GeoObjParser
00047  * geoobj_out
00048  * Envelope
00049  */
00050 
00051 /* Debugging */
00052 #define TRACE(a) elog(NOTICE,a);
00053 
00054 /* For buffer string */
00055 #define MAX_FMTYPE_STRING 20
00056 
00057 /* For the BIG/LITTLE ENDIAN */
00058 #define BYTEORDER 1
00059 
00060 /* For max double string length needs float.h */
00061 #define P_MAXDIG DBL_DIG
00062 #define P_MAXLEN (2*(DBL_DIG+7)+1)
00063 static int  digits8 = DBL_DIG;
00064 
00065 #define max(a,b)        ((a) >  (b) ? (a) : (b))
00066 #define min(a,b)        ((a) <= (b) ? (a) : (b))
00067 #define abs(a)          ((a) <  (0) ? (-a) : (a))
00068 
00069 /* Forward Declarations for things not in geoobj.h */
00070 int getFMTypeCode(char *begin, char **end, int4 *type);
00071 int getFMTypeStr(GeoObj *obj, char *type);
00072 Point *GeoObjPoints(char *begin, char** end, int4 *numpoints);
00073 WKBGeometry *GeoObjParser(char *begin ,char **end, int4 type, int4 *size);
00074 GeoObj *newGeoObj(char *instr);
00075 int makeBox(char *input, char *output);
00076 
00077 /* Input/Output for GeoObj */
00078 
00079 GeoObj *geoobj_in (char *);
00080 char *geoobj_out (GeoObj *);
00081 
00082 /* Return the Bounding Box of a GeoObj as a GeoObj */
00083 GeoObj *envelope (GeoObj *);
00084 
00085 /* 
00086 ** GiST support methods
00087 */
00088 
00089 bool             g_geoobj_consistent(GISTENTRY *entry, GeoObj *query, StrategyNumber strategy);
00090 GISTENTRY *      g_geoobj_compress(GISTENTRY *entry);
00091 GISTENTRY *      g_geoobj_decompress(GISTENTRY *entry);
00092 float *          g_geoobj_penalty(GISTENTRY *origentry, GISTENTRY *newentry, float *result);
00093 GIST_SPLITVEC *  g_geoobj_picksplit(bytea *entryvec, GIST_SPLITVEC *v);
00094 bool             g_geoobj_leaf_consistent(GeoObj *key, GeoObj *query, StrategyNumber strategy);
00095 bool             g_geoobj_internal_consistent(GeoObj *key, GeoObj *query, StrategyNumber strategy);
00096 GeoObj *         g_geoobj_union(bytea *entryvec, int *sizep);
00097 GeoObj *         g_geoobj_binary_union(GeoObj *r1, GeoObj *r2, int *sizep);
00098 bool *           g_geoobj_same(GeoObj *b1, GeoObj *b2, bool *result);
00099 
00100 /*
00101 ** R-tree suport functions
00102 */
00103 bool         geoobj_same(GeoObj *a, GeoObj *b);
00104 bool         geoobj_different(GeoObj *a, GeoObj *b);
00105 bool         geoobj_contains(GeoObj *a, GeoObj *b);
00106 bool         geoobj_contained (GeoObj *a, GeoObj *b);
00107 bool         geoobj_overlap(GeoObj *a, GeoObj *b);
00108 GeoObj *     geoobj_union(GeoObj *a, GeoObj *b);
00109 GeoObj *     geoobj_inter(GeoObj *a, GeoObj *b);
00110 float *      geoobj_size(GeoObj *a);
00111 void         rt_geoobj_size(GeoObj *a, float *sz);
00112 
00113 /*
00114 ** These make no sense for this type, but R-tree wants them
00115 */
00116 bool         geoobj_over_left(GeoObj *a, GeoObj *b);
00117 bool         geoobj_over_right(GeoObj *a, GeoObj *b);
00118 bool         geoobj_left(GeoObj *a, GeoObj *b);
00119 bool         geoobj_right(GeoObj *a, GeoObj *b);
00120 
00121 /* 
00122 ** Auxiliary funxtions
00123 */
00124 static       float distance_1D(float a1, float a2, float b1, float b2);
00125 static       GeoObj *swap_corners (GeoObj *a);
00126 
00136 int getFMTypeCode(char *begin, char **end, int4 *type)
00137 {
00138   int i;
00139   char *cp, *p;
00140   char buffer[MAX_FMTYPE_STRING + 1];
00141 
00142 TRACE("getFMTypeCode")
00143 
00144   p = begin;
00145   *type = 0; 
00146   while(isspace(*p))
00147     p++;
00148   cp = p;
00149   while(isalpha(*p))
00150     p++;
00151   if (p - cp > MAX_FMTYPE_STRING)
00152     return *type;
00153   for(i = 0; cp < p && i < MAX_FMTYPE_STRING; i++, cp++)
00154     buffer[i] = *cp;
00155   buffer[i] = '\0';
00156  
00157   if (strcasecmp(buffer,"POINT") == 0)
00158     *type = FM_POINT;
00159   else if (strcasecmp(buffer,"LINESTRING") == 0)
00160     *type = FM_LINESTRING;
00161   else if (strcasecmp(buffer,"POLYGON") == 0)
00162     *type = FM_POLYGON;
00163   else if (strcasecmp(buffer, "MULTIPOINT") == 0)
00164     *type = FM_MULTIPOINT;
00165   else if (strcasecmp(buffer,"MULTILINESTRING") == 0)
00166     *type = FM_MULTILINESTRING;
00167   else if (strcasecmp(buffer,"MULTIPOLYGON") == 0)
00168     *type = FM_MULTIPOLYGON;
00169   else if (strcasecmp(buffer,"BOX") == 0)
00170     *type = FM_BOX;   
00171     
00172   *end = cp;
00173   return *type;
00174 }
00175 
00176 /* getFMTypeStr: returns the string representation (uppercase ASCII)
00177  * for the GeoObj type, stored as an integer code in the structure.
00178  * On failure it returns zero.  Useful for a function to handle queries
00179  * that want the type string to filter unwanted objects.
00180  */
00181 
00182 int getFMTypeStr(GeoObj *obj, char *type)
00183 {
00184 
00185 TRACE("getFMTypeStr")
00186 
00187   if(obj == NULL || type == NULL)
00188     return 0;
00189 
00190   switch(obj->wkbtype) {
00191   case FM_POINT: 
00192     strncpy(type, "POINT", MAX_FMTYPE_STRING);
00193     return FM_POINT;
00194   case FM_LINESTRING:
00195     strncpy(type, "LINESTRING", MAX_FMTYPE_STRING);
00196     return FM_LINESTRING;
00197   case FM_POLYGON:
00198     strncpy(type,"POLYGON", MAX_FMTYPE_STRING);
00199     return FM_POLYGON;
00200   case FM_MULTILINESTRING:
00201     strncpy(type,"MULTILINESTRING", MAX_FMTYPE_STRING);
00202     return FM_MULTILINESTRING;
00203   case FM_MULTIPOLYGON:
00204     strncpy(type,"MULTIPOLYGON", MAX_FMTYPE_STRING);
00205     return FM_MULTIPOLYGON;
00206   case FM_BOX:
00207     strncpy(type,"BOX", MAX_FMTYPE_STRING);
00208     return FM_BOX;
00209   default:
00210     strncpy(type,"", MAX_FMTYPE_STRING);
00211   }
00212   return 0;
00213 }
00214 
00215 /* GeoObjPoints: retreive an array of points */
00216 
00217 Point *GeoObjPoints(char *begin, char** end, int4 *numpoints)
00218 {
00219   char *p,*cp,*strend;
00220   char buffer[255];
00221   Point *points=NULL;
00222   double x,y,z;
00223 
00224 TRACE("GeoObjPoints")
00225 
00226 
00227   strend = begin;
00228   while(*strend != '\0')
00229     strend++;
00230   
00231   p = begin;
00232   cp = p;  
00233   *numpoints=0;  
00234   
00235 
00236   while(isspace(*p))
00237     p++,cp++;
00238   
00239   if (*p != '(')
00240     elog(ERROR,"Syntax Error: Expecting (");
00241   p++,cp++;
00242 
00243   while(p < strend)
00244   {  
00245     while(isspace(*p))
00246       p++, cp++; 
00247     
00248     x = strtod(p,&cp);
00249     (p != cp) ? p = cp : elog(ERROR, "Numeric conversion error on first coordinate"); 
00250   
00251     while(isspace(*p))
00252       p++, cp++; 
00253          
00254     y = strtod(p,&cp);
00255     (p != cp) ? p = cp : elog(ERROR, "Numeric conversion error on second coordinate");  
00256   
00257     while(isspace(*p))
00258       p++, cp++; 
00259       
00260     z = strtod(p,&cp);
00261     (p != cp) ? p = cp : elog(ERROR, "Numeric conversion error on third coordinate");  
00262     
00263     *numpoints+=1;
00264 
00265     points=realloc(points,sizeof(Point)*(*numpoints));
00266     points[*numpoints-1].x = x;
00267     points[*numpoints-1].y = y;
00268     points[*numpoints-1].z = z;
00269 
00270 sprintf(buffer,"points: %.g %.g %.g",points[*numpoints-1].x,points[*numpoints-1].y,points[*numpoints-1].z);
00271 elog(NOTICE,buffer);
00272 
00273     while(isspace(*p))
00274       p++, cp++;
00275       
00276     if (*p != ',' && *p != ')')
00277       elog(ERROR,"Syntax Error: Expecting ) or ,");
00278     if (*p==')') 
00279     {
00280       p++,cp++;
00281       break;
00282     }
00283     p++,cp++;  
00284 
00285 
00286   }
00287   *end=cp;
00288   return points;
00289   
00290 }
00291 /* GeoObjParser: This function parses an input string and return
00292  * a WKBGeometry that will be stored in a geoobj structure
00293  */
00294  
00295 WKBGeometry *GeoObjParser(char *begin ,char **end, int4 type, int4 *size)
00296 {
00297   WKBGeometry *wkbgeo;
00298   Point *points;
00299   int4 i,numpoints;
00300   char *p,*cp;
00301   char buffer[255];
00302 
00303 TRACE("GeoObjParser")
00304   
00305   *size=0;
00306 
00307   /* initialise the parameters */
00308   switch(type)
00309   {
00310     case FM_POINT:
00311       points = GeoObjPoints(begin,&p,&numpoints);  
00312       if (numpoints!=1) elog(ERROR,"Too many nodes found for type point");
00313       
00314       *size=sizeof(WKBPoint);      
00315       wkbgeo=palloc(*size);
00316       wkbgeo->point.byteOrder = BYTEORDER;  
00317       wkbgeo->point.wkbType = type;
00318       wkbgeo->point.point.x = points[0].x;
00319       wkbgeo->point.point.y = points[0].y;
00320       wkbgeo->point.point.z = points[0].z;
00321       free(points);
00322 
00323       break;
00324       
00325     case FM_LINESTRING:
00326       points = GeoObjPoints(begin,&p,&numpoints);
00327       if (numpoints<=1) elog(ERROR,"Not Enough nodes for type linestring");
00328       *size=sizeof(Point)*(numpoints)+sizeof(WKBLineString);
00329       wkbgeo=palloc(*size);
00330       wkbgeo->linestring.byteOrder = BYTEORDER;
00331       wkbgeo->linestring.wkbType = type;
00332       wkbgeo->linestring.numPoints = numpoints;
00333       for(i=0;i<numpoints;i++)
00334       {
00335         wkbgeo->linestring.points[i].x = points[i].x;
00336         wkbgeo->linestring.points[i].y = points[i].y;
00337         wkbgeo->linestring.points[i].z = points[i].z;         
00338       }   
00339       free(points);
00340       break;
00341 
00342     case FM_POLYGON:
00343       p = begin;
00344       cp = p;
00345       while(isspace(*p))
00346        p++,cp++;
00347   
00348       if (*p != '(')
00349         elog(ERROR,"Syntax Error: Expecting ( after Object type");
00350       p++,cp++;
00351       
00352       points = GeoObjPoints(p,&cp,&numpoints);
00353       (p != cp) ? p = cp : elog(ERROR, "Cannot extract the ring"); 
00354 
00355       while(isspace(*p))
00356        p++,cp++;
00357   
00358       if (*p != ',' && *p != ')')
00359         elog(ERROR,"Syntax Error: Expecting ) or , at the end of the ring");
00360       p++,cp++;
00361 
00362       if (numpoints<=2) elog(ERROR,"Not Enough nodes for type Polygon");
00363       *size=(sizeof(Point)*(numpoints)+sizeof(LinearRing))+sizeof(WKBPolygon);
00364       wkbgeo=palloc(*size);
00365       wkbgeo->polygon.byteOrder = BYTEORDER;
00366       wkbgeo->polygon.wkbType = type;
00367       wkbgeo->polygon.numRings = 1;
00368       wkbgeo->polygon.rings[0].numPoints = numpoints;
00369       
00370       for(i=0;i<numpoints;i++)
00371       {
00372         wkbgeo->polygon.rings[0].points[i].x = points[i].x;
00373         wkbgeo->polygon.rings[0].points[i].y = points[i].y;
00374         wkbgeo->polygon.rings[0].points[i].z = points[i].z;         
00375       }   
00376       free(points);
00377       break;
00378 
00379     case FM_BOX:
00380       points = GeoObjPoints(begin,&p,&numpoints);  
00381       if (numpoints!=2) elog(ERROR,"Two nodes only are required for type box");
00382       
00383       *size=sizeof(WKBBox);      
00384       wkbgeo=palloc(*size);
00385       wkbgeo->box.byteOrder = BYTEORDER; 
00386       wkbgeo->box.wkbType = type;
00387       wkbgeo->box.x[0] = points[0].x;
00388       wkbgeo->box.x[1] = points[0].y;
00389       wkbgeo->box.x[2] = points[0].z;
00390       wkbgeo->box.x[3] = points[1].x;
00391       wkbgeo->box.x[4] = points[1].y;
00392       wkbgeo->box.x[5] = points[1].z;
00393       free(points);
00394 
00395       break;
00396             
00397     default:
00398       elog(ERROR,"Unknown Data Type passed in GeoObjParser");
00399   } 
00400 sprintf(buffer,"size: %d type: %d",*size,type);
00401 elog(NOTICE,buffer); 
00402   *end=p;
00403   return wkbgeo; 
00404   
00405 }
00406 
00407 
00408 /* newGeoObj: This is the main() of functions for creating a new
00409  * GeoObj structure from an input string. It mostly calls other
00410  * functions to do the work. It needs to call a parser to do
00411  * the string copy to remove extra whitespace (TODO).
00412  */
00413 
00414 GeoObj *newGeoObj(char *instr)
00415 {
00416   int4 size, length, type;
00417   char *p;
00418   GeoObj *obj;
00419   WKBGeometry *wkbgeo;
00420   char buffer[255];
00421 
00422 TRACE("newGeoObj")
00423 
00424   /* retreive the type */
00425   p = instr;
00426   if (!getFMTypeCode(instr,&p,&type))
00427     elog(ERROR,"Unknown Data Type!");
00428   /* parse the rest and compact into a WKBGeometry */
00429   instr = p;
00430   wkbgeo = GeoObjParser(instr,&p,type, &length);
00431   
00432   /* store the result into the geoobj */
00433 sprintf(buffer,"length: %d type:%d ",length,type);
00434 elog(NOTICE,buffer);
00435   size = offsetof(GeoObj, data[0]) + length + 1;
00436   obj = palloc(size);
00437   memset(obj,0,size);
00438   obj->size = size;
00439   obj->wkbtype = type;
00440   obj->byte_order = BYTEORDER; 
00441   /* move the wkbgeo structure into the geoobj */
00442   memcpy(obj->data,wkbgeo,length);
00443   
00444   /* free wkbgeo*/
00445   pfree(wkbgeo);
00446   
00447   return obj;
00448 }   
00449 
00450 /*------------------DEPRECIATED----------------------------*/    
00451 /* makeBox: Called with a pointer to an input string
00452  * pointing on or just before the first '(' and a 
00453  * preallocated output string. It returns zero if
00454  * either of the pointers are NULL, otherwise it
00455  * returns the dimensions unless it raises an error.
00456  */
00457 
00458 int makeBox(char *input, char *output)
00459 {
00460   float8 xMin,xMax,yMin,yMax,zMin,zMax,tmp;
00461   int dim;
00462   char *p, *end, *tail;
00463 
00464 TRACE("makeBox")
00465 
00466   if (input == NULL || output == NULL)
00467     return 0;
00468 
00469   p = end = tail = input;
00470   xMin = yMin = zMin = DBL_MAX;
00471   xMax = yMax = zMax = DBL_MIN;
00472   dim = 0;
00473   
00474   while(*end != '\0')
00475     end++;
00476   
00477   while(p < end) {
00478     while(isspace(*p))
00479       p++, tail++;
00480     if(*p != '(')
00481       elog(ERROR, "Unknown character where '(' should have been!");
00482     p++,tail++;
00483     tmp = strtod(p,&tail);
00484     (p != tail) ? p = tail : elog(ERROR, "Numeric conversion error!");
00485     (tmp < xMin) ? xMin = tmp : 1 ;
00486     (tmp > xMax) ? xMax = tmp : 1 ;
00487     while(isspace(*p))
00488       p++, tail++;
00489     if(*p != ',')
00490       elog(ERROR, "Unknown character where ',' should have been!");
00491     p++, tail++;
00492     tmp = strtod(p,&tail);
00493     (p != tail) ? p = tail : elog(ERROR, "Numeric conversion error!");
00494     (tmp < yMin) ? yMin = tmp : 1 ;
00495     (tmp > yMax) ? yMax = tmp : 1 ;
00496     while(isspace(*p))
00497       p++, tail++;
00498     if (dim < 2) { /* checks will catch an error later */
00499       if (*p == ')') {
00500         dim = 2;
00501         zMin = zMax = 0.0;
00502       }
00503       else {
00504         dim = 3;
00505       } 
00506     }
00507     if (dim == 3) {
00508       if(*p != ',')
00509         elog(ERROR, "Unknown character where ',' should have been!");
00510       p++, tail++;
00511       tmp = strtod(p, &tail);
00512       (p != tail) ? p = tail : elog(ERROR, "Numeric conversion error!");
00513       (tmp < zMin) ? zMin = tmp : 1 ;
00514       (tmp > zMax) ? zMax = tmp : 1 ;
00515       while(isspace(*p))
00516         p++, tail++;
00517     }
00518     if(*p != ')')
00519       elog (ERROR, "Unknown character where ')' should have been!");
00520     p++, tail++;
00521     while(isspace(*p))
00522       p++, tail++;
00523   }
00524   sprintf(output,"(%.*g,%.*g,%.*g) (%.*g,%.*g,%.*g)",
00525           digits8, xMin, digits8, yMin, digits8, zMin, 
00526           digits8, xMax, digits8, yMax, digits8, zMax
00527           );
00528 
00529   return dim;
00530 }
00531 
00532 /* geoobj_in: The main entry point for new data from postgresql */
00533 
00534 GeoObj *geoobj_in(char *instr)
00535 {
00536 TRACE("geoobj_in")
00537      return newGeoObj(instr);
00538 }
00539 
00540 
00541 /* geoobj_out: The main output function for the data that postgresql
00542  * calls directly. 
00543  * We have to rebuild the string.
00544  */
00545 
00546 char *geoobj_out (GeoObj *obj)
00547 {
00548   char *outstr,*begin,*coord,*middle;
00549   size_t length;
00550   int4 type;
00551   int4 i;
00552   int4 maxlength;
00553   WKBGeometry *wkbgeo;
00554   char buffer[255];
00555   char *data;
00556 
00557 TRACE("geoobj_out")
00558 
00559   data = obj->data;
00560   wkbgeo = (WKBGeometry *) data;  
00561 
00562 sprintf(buffer,"out: %d ",obj->wkbtype);
00563 elog(NOTICE,buffer);
00564 
00565   switch (obj->wkbtype)
00566   {
00567     case FM_POINT:
00568       maxlength=20*3+MAX_FMTYPE_STRING+4;
00569       outstr=palloc(maxlength);
00570       snprintf(outstr,maxlength,"POINT (%.g %.g %.g)",wkbgeo->point.point.x,wkbgeo->point.point.y,wkbgeo->point.point.z);
00571       break;
00572 
00573     case FM_LINESTRING:
00574       maxlength=wkbgeo->linestring.numPoints*20*3+MAX_FMTYPE_STRING+4;
00575       outstr=palloc(maxlength);
00576       begin=palloc(MAX_FMTYPE_STRING+3);
00577       coord=palloc(20*3);
00578       middle=palloc(maxlength);
00579       
00580       snprintf(begin,MAX_FMTYPE_STRING,"LINESTRING (");
00581       snprintf(middle,maxlength,"");
00582 
00583       for (i=0;i<wkbgeo->linestring.numPoints;i++)
00584       {
00585         snprintf(coord,maxlength,"%.g %.g %.g",wkbgeo->linestring.points[i].x,wkbgeo->linestring.points[i].y,wkbgeo->linestring.points[i].z);
00586         middle=strncat(middle,coord,maxlength);
00587         if (i<wkbgeo->linestring.numPoints-1)
00588         {
00589           middle=strncat(middle,", ",maxlength);
00590         } 
00591       }
00592       snprintf(outstr,maxlength,"%s%s)",begin,middle);
00593       pfree(begin);
00594       pfree(coord);
00595       pfree(middle);
00596       break;
00597       
00598     case FM_POLYGON:
00599       maxlength=wkbgeo->polygon.numRings*wkbgeo->polygon.rings[0].numPoints*20*3+MAX_FMTYPE_STRING+4;
00600       outstr=palloc(maxlength);
00601       begin=palloc(MAX_FMTYPE_STRING+3);
00602       coord=palloc(20*3);
00603       middle=palloc(maxlength);
00604       
00605       snprintf(begin,MAX_FMTYPE_STRING,"POLYGON ((");
00606       snprintf(middle,maxlength,"");
00607 
00608       for (i=0;i<wkbgeo->polygon.rings[0].numPoints;i++)
00609       {
00610         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);
00611         middle=strncat(middle,coord,maxlength);
00612         if (i<wkbgeo->polygon.rings[0].numPoints-1)
00613         {
00614           middle=strncat(middle,", ",maxlength);
00615         } 
00616       }
00617       snprintf(outstr,maxlength,"%s%s))",begin,middle);
00618       pfree(begin);
00619       pfree(coord);
00620       pfree(middle);
00621       break;
00622 
00623     case FM_BOX:
00624       maxlength=20*6+MAX_FMTYPE_STRING+4;
00625       outstr=palloc(maxlength);
00626       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]);
00627       break;
00628             
00629     default: 
00630       elog(ERROR,"Unknown type in geoobj_out");
00631   }
00632 elog(NOTICE,outstr);  
00633   return outstr;
00634 }
00635 
00636 /* ----------------------DEPRECIATED --------------------------*/
00637 /* box_out: Makes a box for any object as a new GeoObj, but
00638  * handles the case where the type is already a box specially.
00639  */
00640 
00641 GeoObj *bounds (GeoObj *obj)
00642 {
00643   char    *outstr, *cp;
00644   int     dims, length;
00645 
00646 TRACE("bounds")
00647 
00648   if (obj == NULL)
00649     elog (ERROR,"box_out was passed a NULL object!");
00650   
00651   if (obj->wkbtype == FM_BOX)
00652     return obj;
00653   else
00654     dims = 3;
00655 
00656   length = dims * P_MAXLEN * 2 + sizeof("BOX ");
00657   outstr = palloc(length + 1);
00658   memset(outstr,0,length + 1);
00659   strcpy(outstr,"BOX ");
00660   cp = outstr;
00661   while(*cp != '\0')
00662     cp++;
00663   if (!makeBox(strchr(obj->data,'('),cp))
00664     elog (ERROR, "makeBox failed!\n");
00665 
00666   return newGeoObj(outstr);
00667 }
00668 
00669 /* Envelope: Returns the Bounding Box of any geoobj as a BOX type
00670  *
00671  */
00672 
00673 GeoObj *envelope (GeoObj *obj)
00674 {
00675   GeoObj *result;
00676   char strbox[20*6+MAX_FMTYPE_STRING+4];
00677   double minx;
00678   double miny;
00679   double minz;
00680   double maxx;
00681   double maxy;
00682   double maxz;
00683   int4 i;
00684   char *data;
00685   WKBGeometry *wkbgeo;
00686   char buffer[255];
00687 
00688 TRACE("envelope")
00689   
00690 sprintf(buffer,"envelope: %d size:%d",obj->wkbtype,obj->size);  
00691 elog(NOTICE,buffer);
00692   data = obj->data;
00693   wkbgeo = (WKBGeometry *) data;
00694   
00695   switch (obj->wkbtype)
00696   {
00697     case FM_POINT:
00698 elog(NOTICE,"point");
00699       minx=maxx=wkbgeo->point.point.x;
00700       miny=maxy=wkbgeo->point.point.y;
00701       minz=maxz=wkbgeo->point.point.z;
00702       break;
00703 
00704     case FM_LINESTRING:
00705       minx=maxx=wkbgeo->linestring.points[0].x;
00706       miny=maxy=wkbgeo->linestring.points[0].y;
00707       minz=maxz=wkbgeo->linestring.points[0].z;
00708       
00709       for (i=0;i<wkbgeo->linestring.numPoints;i++)
00710       {
00711         if (wkbgeo->linestring.points[i].x<minx) minx=wkbgeo->linestring.points[i].x;
00712         if (wkbgeo->linestring.points[i].y<miny) miny=wkbgeo->linestring.points[i].y;
00713         if (wkbgeo->linestring.points[i].z<minz) minz=wkbgeo->linestring.points[i].z;
00714         if (wkbgeo->linestring.points[i].x>maxx) maxx=wkbgeo->linestring.points[i].x;
00715         if (wkbgeo->linestring.points[i].y>maxy) maxy=wkbgeo->linestring.points[i].y;
00716         if (wkbgeo->linestring.points[i].z>maxz) maxz=wkbgeo->linestring.points[i].z;
00717       }      
00718       break;
00719       
00720     case FM_POLYGON:
00721       minx=maxx=wkbgeo->polygon.rings[0].points[0].x;
00722       miny=maxy=wkbgeo->polygon.rings[0].points[0].y;
00723       minz=maxz=wkbgeo->polygon.rings[0].points[0].z;
00724       
00725       for (i=0;i<wkbgeo->polygon.rings[0].numPoints;i++)
00726       {
00727         if (wkbgeo->polygon.rings[0].points[i].x<minx) minx=wkbgeo->polygon.rings[0].points[i].x;
00728         if (wkbgeo->polygon.rings[0].points[i].y<miny) miny=wkbgeo->polygon.rings[0].points[i].y;
00729         if (wkbgeo->polygon.rings[0].points[i].z<minz) minz=wkbgeo->polygon.rings[0].points[i].z;
00730         if (wkbgeo->polygon.rings[0].points[i].x>maxx) maxx=wkbgeo->polygon.rings[0].points[i].x;
00731         if (wkbgeo->polygon.rings[0].points[i].y>maxy) maxy=wkbgeo->polygon.rings[0].points[i].y;
00732         if (wkbgeo->polygon.rings[0].points[i].z>maxz) maxz=wkbgeo->polygon.rings[0].points[i].z;
00733       }
00734       break;  
00735       
00736     case FM_BOX:    
00737       minx=wkbgeo->box.x[0];
00738       miny=wkbgeo->box.x[1];
00739       minz=wkbgeo->box.x[2];
00740       maxx=wkbgeo->box.x[3];
00741       maxy=wkbgeo->box.x[4];
00742       maxz=wkbgeo->box.x[5];
00743       break;
00744             
00745     default:
00746       elog(ERROR,"Unknown type in envelope");
00747   }
00748   sprintf(strbox,"BOX (%.g %.g %.g, %.g %.g %.g)",minx,miny,minz,maxx,maxy,maxz);
00749 elog(NOTICE,strbox);
00750   result=newGeoObj(strbox);
00751   return(result);
00752 }
00753 /*****************************************************************************
00754  *                         GiST functions
00755  *****************************************************************************/
00756 /*
00757 ** The GiST Consistent method for boxes
00758 ** Should return false if for all data items x below entry,
00759 ** the predicate x op query == FALSE, where op is the oper
00760 ** corresponding to strategy in the pg_amop table.
00761 */
00762 
00763 bool g_geoobj_consistent(GISTENTRY *entry, GeoObj *query, StrategyNumber strategy)
00764 {
00765 TRACE("g_geoobj_consistent")
00766 
00767     /*
00768     ** if entry is not leaf, use g_geobj_internal_consistent,
00769     ** else use g_geoobj_leaf_consistent
00770     */
00771     if (GIST_LEAF(entry))
00772       return(g_geoobj_leaf_consistent((GeoObj *)(entry->pred), query, strategy));
00773     else
00774       return(g_geoobj_internal_consistent((GeoObj *)(entry->pred), query, strategy));
00775 }
00776 
00777 /*
00778 ** GiST Compress and Decompress methods for geoobj
00779 ** Compress object into its bounding box
00780 ** decompres the bounding box into itself
00781 */
00782 
00783 GISTENTRY *g_geoobj_compress(GISTENTRY *entry)
00784 {
00785   char buffer[255];
00786   GeoObj *obj;
00787 TRACE("g_geoobj_compress")
00788 obj=(GeoObj *)entry->pred;
00789 sprintf(buffer,"type: %d",obj->wkbtype);
00790 elog(NOTICE,buffer);
00791 elog(NOTICE,geoobj_out(obj));
00792     entry->pred=(char *)envelope((GeoObj *)entry->pred);
00793 
00794     return(entry);
00795 }
00796 
00797 GISTENTRY *g_geoobj_decompress(GISTENTRY *entry)
00798 {
00799   char buffer[255];
00800   GeoObj *obj;
00801 TRACE("g_geoobj_decompress")
00802 obj=(GeoObj *)entry->pred;
00803 sprintf(buffer,"type: %d",obj->wkbtype);
00804 elog(NOTICE,buffer);
00805 elog(NOTICE,geoobj_out(obj));
00806     return(entry);
00807 }
00808 
00809 /* 
00810 ** SUPPORT ROUTINES
00811 */
00812 bool g_geoobj_leaf_consistent(GeoObj *key, GeoObj *query, StrategyNumber strategy)
00813 {
00814   bool retval;
00815 
00816 TRACE("g_geoobj_leaf_consistent")
00817   /*
00818   fprintf(stderr, "leaf_consistent, %d\n", strategy);
00819   */
00820   switch(strategy) {
00821   case RTLeftStrategyNumber:
00822     retval = (bool)geoobj_left(key, query);
00823     break;
00824   case RTOverLeftStrategyNumber:
00825     retval = (bool)geoobj_over_left(key,query);
00826     break;
00827   case RTOverlapStrategyNumber:
00828     retval = (bool)geoobj_overlap(key, query);
00829     break;
00830   case RTOverRightStrategyNumber:
00831     retval = (bool)geoobj_over_right(key, query);
00832     break;
00833   case RTRightStrategyNumber:
00834     retval = (bool)geoobj_right(key, query);
00835     break;
00836   case RTSameStrategyNumber:
00837     retval = (bool)geoobj_same(key, query);
00838     break;
00839   case RTContainsStrategyNumber:
00840     retval = (bool)geoobj_contains(key, query);
00841     break;
00842   case RTContainedByStrategyNumber:
00843     retval = (bool)geoobj_contained(key,query);
00844     break;
00845   default:
00846     retval = FALSE;
00847   }
00848   return(retval);
00849 }
00850 
00851 bool g_geoobj_internal_consistent(GeoObj *key, GeoObj *query, StrategyNumber strategy)
00852 {
00853   bool retval;
00854 
00855 TRACE("g_geoobj_internal_consistent")
00856   
00857   /*
00858   fprintf(stderr, "internal_consistent, %d\n", strategy);
00859   */
00860   switch(strategy) {
00861   case RTLeftStrategyNumber:
00862   case RTOverLeftStrategyNumber:
00863     retval = (bool)geoobj_over_left(key,query);
00864     break;
00865   case RTOverlapStrategyNumber:
00866     retval = (bool)geoobj_overlap(key, query);
00867     break;
00868   case RTOverRightStrategyNumber:
00869   case RTRightStrategyNumber:
00870     retval = (bool)geoobj_right(key, query);
00871     break;
00872   case RTSameStrategyNumber:
00873   case RTContainsStrategyNumber:
00874     retval = (bool)geoobj_contains(key, query);
00875     break;
00876   case RTContainedByStrategyNumber:
00877     retval = (bool)geoobj_overlap(key, query);
00878     break;
00879   default:
00880     retval = FALSE;
00881     }
00882   return(retval);
00883 }
00884 /*
00885 ** End SUPPORT ROUTINES
00886 */
00887 
00888 
00889 /*
00890 ** The GiST Union method for boxes
00891 ** returns the minimal bounding box that encloses all the entries in entryvec
00892 */
00893 GeoObj *g_geoobj_union(bytea *entryvec, int *sizep)
00894 {
00895     int numranges, i;
00896     GeoObj *out = (GeoObj *)NULL;
00897     GeoObj *tmp;
00898 TRACE("g_geoobj_union")
00899     /*
00900     fprintf(stderr, "union\n");
00901     */
00902     numranges = (VARSIZE(entryvec) - VARHDRSZ)/sizeof(GISTENTRY); 
00903     tmp = (GeoObj *)(((GISTENTRY *)(VARDATA(entryvec)))[0]).pred;
00904     /*
00905     *sizep = sizeof(NDBOX); -- NDBOX has variable size
00906     */
00907     *sizep = tmp->size;
00908 
00909     for (i = 1; i < numranges; i++) {
00910       out = g_geoobj_binary_union(tmp, (GeoObj *)
00911                                  (((GISTENTRY *)(VARDATA(entryvec)))[i]).pred,
00912                                  sizep);
00913       /*
00914         fprintf(stderr, "\t%s ^ %s -> %s\n", cube_out(tmp), cube_out((NDBOX *)(((GISTENTRY *)(VARDATA(entryvec)))[i]).pred), cube_out(out));
00915       */
00916       if (i > 1) pfree(tmp);
00917       tmp = out;
00918     }
00919 
00920     return(out);
00921 }
00922 
00923 /*
00924 ** The GiST Penalty method for boxes
00925 ** As in the R-tree paper, we use change in area as our penalty metric
00926 */
00927 float *g_geoobj_penalty(GISTENTRY *origentry, GISTENTRY *newentry, float *result)
00928 {
00929     Datum ud;
00930     float tmp1, tmp2;
00931 
00932 TRACE("g_geoobj_penalty")
00933     
00934     ud = (Datum)geoobj_union((GeoObj *)(origentry->pred), (GeoObj *)(newentry->pred));
00935     rt_geoobj_size((GeoObj *)ud, &tmp1);
00936     rt_geoobj_size((GeoObj *)(origentry->pred), &tmp2);
00937     *result = tmp1 - tmp2;
00938     pfree((char *)ud);
00939     /*
00940     fprintf(stderr, "penalty\n");
00941     fprintf(stderr, "\t%g\n", *result);
00942     */
00943     return(result);
00944 }
00945 
00946 
00947 
00948 /*
00949 ** The GiST PickSplit method for boxes
00950 ** We use Guttman's poly time split algorithm 
00951 */
00952 GIST_SPLITVEC *g_geoobj_picksplit(bytea *entryvec, GIST_SPLITVEC *v)
00953 {
00954     OffsetNumber i, j;
00955     GeoObj *datum_alpha, *datum_beta;
00956     GeoObj *datum_l, *datum_r;
00957     GeoObj *union_d, *union_dl, *union_dr;
00958     GeoObj *inter_d;
00959     bool firsttime;
00960     float size_alpha, size_beta, size_union, size_inter;
00961     float size_waste, waste;
00962     float size_l, size_r;
00963     int nbytes;
00964     OffsetNumber seed_1 = 0, seed_2 = 0;
00965     OffsetNumber *left, *right;
00966     OffsetNumber maxoff;
00967 
00968 TRACE("g_geoobj_picksplit")
00969 
00970     /*
00971     fprintf(stderr, "picksplit\n");
00972     */
00973     maxoff = ((VARSIZE(entryvec) - VARHDRSZ)/sizeof(GISTENTRY)) - 2;
00974     nbytes =  (maxoff + 2) * sizeof(OffsetNumber);
00975     v->spl_left = (OffsetNumber *) palloc(nbytes);
00976     v->spl_right = (OffsetNumber *) palloc(nbytes);
00977     
00978     firsttime = true;
00979     waste = 0.0;
00980     
00981     for (i = FirstOffsetNumber; i < maxoff; i = OffsetNumberNext(i)) 
00982     {
00983             datum_alpha = (GeoObj *)(((GISTENTRY *)(VARDATA(entryvec)))[i].pred);
00984             for (j = OffsetNumberNext(i); j <= maxoff; j = OffsetNumberNext(j)) 
00985             {
00986               datum_beta = (GeoObj *)(((GISTENTRY *)(VARDATA(entryvec)))[j].pred);
00987             
00988               /* compute the wasted space by unioning these guys */
00989             /* size_waste = size_union - size_inter; */
00990             union_d = (GeoObj *)geoobj_union(datum_alpha, datum_beta);
00991             rt_geoobj_size(union_d, &size_union);
00992             inter_d = (GeoObj *)geoobj_inter(datum_alpha, datum_beta);
00993             rt_geoobj_size(inter_d, &size_inter);
00994             size_waste = size_union - size_inter;
00995             
00996             pfree(union_d);
00997             
00998             if (inter_d != (GeoObj *) NULL)
00999                 pfree(inter_d);
01000             
01001               /*
01002              *  are these a more promising split than what we've
01003              *  already seen?
01004              */
01005             
01006               if (size_waste > waste || firsttime) 
01007               {
01008                 waste = size_waste;
01009           seed_1 = i;
01010                       seed_2 = j;
01011                       firsttime = false;
01012               }
01013             }
01014     }
01015     
01016     left = v->spl_left;
01017     v->spl_nleft = 0;
01018     right = v->spl_right;
01019     v->spl_nright = 0;
01020     
01021     datum_alpha = (GeoObj *)(((GISTENTRY *)(VARDATA(entryvec)))[seed_1].pred);
01022     datum_l = (GeoObj *)geoobj_union(datum_alpha, datum_alpha);
01023     rt_geoobj_size((GeoObj *)datum_l, &size_l);
01024     datum_beta = (GeoObj *)(((GISTENTRY *)(VARDATA(entryvec)))[seed_2].pred);;
01025     datum_r = (GeoObj *)geoobj_union(datum_beta, datum_beta);
01026     rt_geoobj_size((GeoObj *)datum_r, &size_r);
01027     
01028     /*
01029      *  Now split up the regions between the two seeds.  An important
01030      *  property of this split algorithm is that the split vector v
01031      *  has the indices of items to be split in order in its left and
01032      *  right vectors.  We exploit this property by doing a merge in
01033      *  the code that actually splits the page.
01034      *
01035      *  For efficiency, we also place the new index tuple in this loop.
01036      *  This is handled at the very end, when we have placed all the
01037      *  existing tuples and i == maxoff + 1.
01038      */
01039     
01040     maxoff = OffsetNumberNext(maxoff);
01041     for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i)) 
01042     {
01043         
01044             /*
01045          *  If we've already decided where to place this item, just
01046          *  put it on the right list.  Otherwise, we need to figure
01047          *  out which page needs the least enlargement in order to
01048          *  store the item.
01049          */
01050         
01051             if (i == seed_1) 
01052             {
01053               *left++ = i;
01054               v->spl_nleft++;
01055               continue;
01056             } else if (i == seed_2) 
01057             {
01058               *right++ = i;
01059               v->spl_nright++;
01060               continue;
01061             }
01062         
01063             /* okay, which page needs least enlargement? */ 
01064         datum_alpha = (GeoObj *)(((GISTENTRY *)(VARDATA(entryvec)))[i].pred);
01065         union_dl = (GeoObj *)geoobj_union(datum_l, datum_alpha);
01066         union_dr = (GeoObj *)geoobj_union(datum_r, datum_alpha);
01067         rt_geoobj_size((GeoObj *)union_dl, &size_alpha);
01068         rt_geoobj_size((GeoObj *)union_dr, &size_beta);
01069         
01070         /* pick which page to add it to */
01071         if (size_alpha - size_l < size_beta - size_r) 
01072         {
01073               pfree(datum_l);
01074             pfree(union_dr);
01075             datum_l = union_dl;
01076             size_l = size_alpha;
01077             *left++ = i;
01078             v->spl_nleft++;
01079         } else 
01080         {
01081               pfree(datum_r);
01082             pfree(union_dl);
01083             datum_r = union_dr;
01084             size_r = size_alpha;
01085             *right++ = i;
01086             v->spl_nright++;
01087             }
01088     }
01089     *left = *right = FirstOffsetNumber; /* sentinel value, see dosplit() */
01090     
01091     v->spl_ldatum = (char *)datum_l;
01092     v->spl_rdatum = (char *)datum_r;
01093 
01094     return v;
01095 }
01096 
01097 /*
01098 ** Equality method
01099 */
01100 bool *g_geoobj_same(GeoObj *b1, GeoObj *b2, bool *result)
01101 {
01102   TRACE("g_geoobj_same")
01103   if (geoobj_same(b1, b2))
01104     *result = TRUE;
01105   else *result = FALSE;
01106   /*
01107   fprintf(stderr, "same: %s\n", (*result ? "TRUE" : "FALSE" ));
01108   */
01109   return(result);
01110 }
01111 
01112 GeoObj *g_geoobj_binary_union(GeoObj *r1, GeoObj *r2, int *sizep)
01113 {
01114     GeoObj *retval;
01115 
01116 TRACE("g_geoobj_binary_union")
01117     retval = geoobj_union(r1, r2);
01118     *sizep = retval->size;
01119 
01120     return (retval);
01121 }
01122 
01123 
01124 /* cube_union */
01125 GeoObj *geoobj_union(GeoObj *box_a, GeoObj *box_b)
01126 {
01127   int i;
01128   GeoObj *result;
01129   GeoObj *a = swap_corners(box_a);
01130   GeoObj *b = swap_corners(box_b);
01131   char *data;
01132   WKBGeometry *wkbgeoresult;
01133   WKBGeometry *wkbgeoa;
01134   WKBGeometry *wkbgeob;
01135 
01136 TRACE("geoobj_union")  
01137   data=a->data;
01138   wkbgeoa=(WKBGeometry *)data;
01139 
01140   data=b->data;
01141   wkbgeob=(WKBGeometry *)data;
01142   
01143   result = palloc(a->size);
01144   result->size = a->size;
01145   result->byte_order = BYTEORDER;
01146   result->wkbtype = FM_BOX;
01147   data=result->data;
01148   wkbgeoresult=(WKBGeometry *)data;
01149   wkbgeoresult->box.byteOrder = BYTEORDER;
01150   wkbgeoresult->box.wkbType = FM_BOX;
01151 
01152   
01153   /* use the potentially smaller of the two boxes (b) to fill in 
01154      the result, padding absent dimensions with zeroes*/
01155   for ( i = 0; i < 3; i++ ) {
01156     wkbgeoresult->box.x[i] = wkbgeob->box.x[i];
01157     wkbgeoresult->box.x[i + 3] = wkbgeob->box.x[i + 3];
01158   }
01159 
01160     
01161   /* compute the union */
01162   for ( i = 0; i < 3; i++ ) {
01163     wkbgeoresult->box.x[i] = min(wkbgeoa->box.x[i], wkbgeoresult->box.x[i]);
01164   }
01165   for ( i = 3; i < 6; i++ ) {
01166     wkbgeoresult->box.x[i] = max(wkbgeoa->box.x[i], wkbgeoresult->box.x[i]);
01167   }
01168 
01169   pfree(a);
01170   pfree(b);
01171 
01172   return(result);
01173 }
01174 
01175 /* cube_inter */
01176 GeoObj *geoobj_inter(GeoObj *box_a, GeoObj *box_b)
01177 {
01178   int i;
01179   GeoObj * result;
01180   GeoObj *a = swap_corners(box_a);
01181   GeoObj *b = swap_corners(box_b);
01182   char *data;
01183   WKBGeometry *wkbgeoresult;
01184   WKBGeometry *wkbgeoa;
01185   WKBGeometry *wkbgeob;
01186 
01187 TRACE("geoobj_inter")  
01188   data=a->data;
01189   wkbgeoa=(WKBGeometry *)data;
01190 
01191   data=b->data;
01192   wkbgeob=(WKBGeometry *)data;
01193     
01194   result = palloc(a->size);
01195   result->size = a->size;
01196   result->byte_order = BYTEORDER;
01197   result->wkbtype = FM_BOX;
01198   data=result->data;
01199   wkbgeoresult=(WKBGeometry *)data;
01200   wkbgeoresult->box.byteOrder = BYTEORDER;
01201   wkbgeoresult->box.wkbType = FM_BOX;
01202 
01203   /* use the potentially  smaller of the two boxes (b) to fill in 
01204      the result, padding absent dimensions with zeroes*/
01205   for ( i = 0; i < 3; i++ ) {
01206     wkbgeoresult->box.x[i] = wkbgeob->box.x[i];
01207     wkbgeoresult->box.x[i + 3] = wkbgeob->box.x[i + 3];
01208   }
01209     
01210   /* compute the intersection */
01211   for ( i = 0; i < 3; i++ ) {
01212     wkbgeoresult->box.x[i] = max(wkbgeoa->box.x[i], wkbgeoresult->box.x[i]);
01213   }
01214   for ( i = 3; i < 6; i++ ) {
01215     wkbgeoresult->box.x[i] = min(wkbgeoa->box.x[i], wkbgeoresult->box.x[i]);
01216   }
01217   
01218   pfree(a);
01219   pfree(b);
01220 
01221   /* Is it OK to return a non-null intersection for non-overlapping boxes? */
01222   return(result);
01223 }
01224 
01225 /* cube_size */
01226 float *geoobj_size(GeoObj *a)
01227 {
01228   int i,j;
01229   float *result;
01230   char *data;
01231   WKBGeometry *wkbgeoa;
01232 
01233 TRACE("geoobj_size")  
01234   data=a->data;
01235   wkbgeoa=(WKBGeometry *)data;
01236   
01237   result = (float *) palloc(sizeof(float));
01238   
01239   *result = 1.0;
01240   for ( i = 0, j = 3; i < 3; i++,j++ ) {
01241     *result=(*result)*abs((wkbgeoa->box.x[j] - wkbgeoa->box.x[i]));
01242   }
01243   
01244   return(result);
01245 }
01246 
01247 void rt_geoobj_size(GeoObj *a, float *size)
01248 {
01249   int i,j;
01250   char *data;
01251   WKBGeometry *wkbgeoa;
01252 
01253 TRACE("rt_geoobj_size")  
01254   data=a->data;
01255   wkbgeoa=(WKBGeometry *)data;
01256   
01257   if (a == (GeoObj *) NULL)
01258     *size = 0.0;
01259   else {
01260     *size = 1.0;
01261     for ( i = 0, j = 3; i < 3; i++,j++ ) {
01262       *size=(*size)*abs((wkbgeoa->box.x[j] - wkbgeoa->box.x[i]));
01263     }
01264   }
01265   return;
01266 }
01267 
01268 /* The following four methods compare the projections of the boxes
01269    onto the 0-th coordinate axis. These methods are useless for dimensions
01270    larger than 2, but it seems that R-tree requires all its strategies
01271    map to real functions that return something */
01272 
01273 /*  is the right edge of (a) located to the left of 
01274     the right edge of (b)? */
01275 bool geoobj_over_left(GeoObj *box_a, GeoObj *box_b)
01276 {
01277   GeoObj *a;
01278   GeoObj *b;
01279   char *data;
01280   WKBGeometry *wkbgeoa;
01281   WKBGeometry *wkbgeob;
01282 
01283 TRACE("geoobj_over_left")  
01284   if ( (box_a==NULL) || (box_b==NULL) )
01285     return(FALSE);
01286 
01287   a = swap_corners(box_a);
01288   b = swap_corners(box_b);
01289 
01290   data=a->data;
01291   wkbgeoa=(WKBGeometry *)data;
01292 
01293   data=b->data;
01294   wkbgeob=(WKBGeometry *)data;
01295   
01296   return( wkbgeoa->box.x[3 - 1] <= wkbgeob->box.x[3 - 1] && !geoobj_left(a, b) && !geoobj_right(a, b) );
01297 }
01298 
01299 /*  is the left edge of (a) located to the right of 
01300     the left edge of (b)? */
01301 bool geoobj_over_right(GeoObj *box_a, GeoObj *box_b)
01302 {
01303   GeoObj *a;
01304   GeoObj *b;
01305   char *data;
01306   WKBGeometry *wkbgeoa;
01307   WKBGeometry *wkbgeob;
01308 
01309 TRACE("geoobj_over_right")  
01310   if ( (box_a==NULL) || (box_b==NULL) )
01311     return(FALSE);
01312 
01313   a = swap_corners(box_a);
01314   b = swap_corners(box_b);
01315 
01316   data=a->data;
01317   wkbgeoa=(WKBGeometry *)data;
01318 
01319   data=b->data;
01320   wkbgeob=(WKBGeometry *)data; 
01321 
01322   return( wkbgeoa->box.x[3 - 1] >= wkbgeob->box.x[3 - 1] && !geoobj_left(a, b) && !geoobj_right(a, b) );
01323 }
01324 
01325 
01326 /* return 'true' if the projection of 'a' is
01327    entirely on the left of the projection of 'b' */
01328 bool geoobj_left(GeoObj *box_a, GeoObj *box_b)
01329 {
01330   GeoObj *a;
01331   GeoObj *b;
01332   char *data;
01333   WKBGeometry *wkbgeoa;
01334   WKBGeometry *wkbgeob;
01335 
01336 TRACE("geoobj_left")    
01337   if ( (box_a==NULL) || (box_b==NULL) )
01338     return(FALSE);
01339 
01340   a = swap_corners(box_a);
01341   b = swap_corners(box_b);
01342   
01343   data=a->data;
01344   wkbgeoa=(WKBGeometry *)data;
01345 
01346   data=b->data;
01347   wkbgeob=(WKBGeometry *)data;
01348   
01349   return( wkbgeoa->box.x[3 - 1] < wkbgeob->box.x[0]);
01350 }
01351 
01352 /* return 'true' if the projection of 'a' is
01353    entirely on the right  of the projection of 'b' */
01354 bool geoobj_right(GeoObj *box_a, GeoObj *box_b)
01355 {
01356   GeoObj *a;
01357   GeoObj *b;
01358   char *data;
01359   WKBGeometry *wkbgeoa;
01360   WKBGeometry *wkbgeob;
01361 
01362 TRACE("geoobj_right")    
01363   if ( (box_a==NULL) || (box_b==NULL) )
01364     return(FALSE);
01365 
01366   a = swap_corners(box_a);
01367   b = swap_corners(box_b);
01368 
01369   data=a->data;
01370   wkbgeoa=(WKBGeometry *)data;
01371 
01372   data=b->data;
01373   wkbgeob=(WKBGeometry *)data;
01374 
01375   return( wkbgeoa->box.x[0] > wkbgeob->box.x[3 - 1]);
01376 }
01377 
01378 /* make up a metric in which one box will be 'lower' than the other
01379    -- this can be useful for srting and to determine uniqueness */
01380 bool geoobj_lt(GeoObj *box_a, GeoObj *box_b)
01381 {
01382   int i;
01383   GeoObj *a;
01384   GeoObj *b;
01385   char *data;
01386   WKBGeometry *wkbgeoa;
01387   WKBGeometry *wkbgeob;
01388 
01389 TRACE("geoobj_lt")    
01390   if ( (box_a==NULL) || (box_b==NULL) )
01391     return(FALSE);
01392 
01393   a = swap_corners(box_a);
01394   b = swap_corners(box_b);
01395 
01396   data=a->data;
01397   wkbgeoa=(WKBGeometry *)data;
01398 
01399   data=b->data;
01400   wkbgeob=(WKBGeometry *)data;
01401 
01402   /* if all common dimensions are equal, the cube with more dimensions wins */
01403   if ( geoobj_same(a, b) ) {
01404     return(FALSE);
01405   }
01406 
01407   /* compare the common dimensions */
01408   for ( i = 0; i < 3; i++ ) {
01409     if ( wkbgeoa->box.x[i] > wkbgeob->box.x[i] )
01410       return(FALSE); 
01411     if ( wkbgeoa->box.x[i] < wkbgeob->box.x[i] )
01412       return(TRUE); 
01413   }
01414   for ( i = 0; i < 3; i++ ) {
01415     if ( wkbgeoa->box.x[i + 3] > wkbgeob->box.x[i + 3] )
01416       return(FALSE); 
01417     if ( wkbgeoa->box.x[i +3] < wkbgeob->box.x[i + 3] )
01418       return(TRUE); 
01419   }
01420 
01421   return(FALSE);
01422 }
01423 
01424 
01425 bool geoobj_gt(GeoObj *box_a, GeoObj *box_b)
01426 {
01427   int i;
01428   GeoObj *a;
01429   GeoObj *b;
01430   char *data;
01431   WKBGeometry *wkbgeoa;
01432   WKBGeometry *wkbgeob;
01433 
01434 TRACE("geoobj_gt")    
01435   if ( (box_a==NULL) || (box_b==NULL) )
01436     return(FALSE);
01437 
01438   a = swap_corners(box_a);
01439   b = swap_corners(box_b);
01440 
01441   data=a->data;
01442   wkbgeoa=(WKBGeometry *)data;
01443 
01444   data=b->data;
01445   wkbgeob=(WKBGeometry *)data;
01446   
01447   /* if all common dimensions are equal, the cube with more dimensions wins */
01448   if ( geoobj_same(a, b) ) 
01449   {
01450     return(FALSE);
01451   }
01452 
01453   /* compare the common dimensions */
01454   for ( i = 0; i < 3; i++ ) 
01455   {
01456     if ( wkbgeoa->box.x[i] < wkbgeob->box.x[i] )
01457       return(FALSE); 
01458     if ( wkbgeoa->box.x[i] > wkbgeob->box.x[i] )
01459       return(TRUE); 
01460   }
01461   for ( i = 0; i < 3; i++ ) 
01462   {
01463     if ( wkbgeoa->box.x[i + 3] < wkbgeob->box.x[i + 3] )
01464       return(FALSE); 
01465     if ( wkbgeoa->box.x[i + 3] > wkbgeob->box.x[i + 3] )
01466       return(TRUE); 
01467   }
01468 
01469   return(FALSE);
01470 }
01471 
01472 
01473 /* Equal */
01474 bool geoobj_same(GeoObj *box_a, GeoObj *box_b)
01475 {
01476   int i;
01477   GeoObj *a;
01478   GeoObj *b;
01479   char *data;
01480   WKBGeometry *wkbgeoa;
01481   WKBGeometry *wkbgeob;
01482 
01483 TRACE("geoobj_same")    
01484   if ( (box_a==NULL) || (box_b==NULL) )
01485     return(FALSE);
01486 
01487   a = swap_corners(box_a);
01488   b = swap_corners(box_b);
01489   
01490   data=a->data;
01491   wkbgeoa=(WKBGeometry *)data;
01492 
01493   data=b->data;
01494   wkbgeob=(WKBGeometry *)data;
01495 
01496   for ( i = 0; i < 3; i++ ) {
01497     if ( wkbgeoa->box.x[i] != wkbgeob->box.x[i] )
01498       return(FALSE);
01499     if ( wkbgeoa->box.x[i + 3] != wkbgeob->box.x[i + 3] )
01500       return(FALSE);
01501   }
01502 
01503   pfree(a);
01504   pfree(b);
01505 
01506   return(TRUE);
01507 }
01508 
01509 /* Different */
01510 bool geoobj_different(GeoObj *box_a, GeoObj *box_b)
01511 {
01512 TRACE("geoobj_different")
01513   return(!geoobj_same(box_a, box_b));
01514 }
01515 
01516 
01517 /* Contains */
01518 /* Box(A) CONTAINS Box(B) IFF pt(A) < pt(B) */
01519 bool geoobj_contains(GeoObj *box_a, GeoObj *box_b)
01520 {
01521   int i;
01522   GeoObj *a;
01523   GeoObj *b;
01524   char *data;
01525   WKBGeometry *wkbgeoa;
01526   WKBGeometry *wkbgeob;
01527 
01528 TRACE("geoobj_contains")  
01529   if ( (box_a==NULL) || (box_b==NULL) )
01530     return(FALSE);
01531 
01532   a = swap_corners(box_a);
01533   b = swap_corners(box_b);
01534 
01535   data=a->data;
01536   wkbgeoa=(WKBGeometry *)data;
01537 
01538   data=b->data;
01539   wkbgeob=(WKBGeometry *)data;
01540   
01541   /* Can't care less about the excess dimensions of (a), if any */
01542   for ( i = 0; i < 3; i++ ) {
01543     if ( wkbgeoa->box.x[i] > wkbgeob->box.x[i] )
01544       return(FALSE);
01545     if ( wkbgeoa->box.x[i + 3] < wkbgeob->box.x[i + 3] )
01546       return(FALSE);
01547   }
01548 
01549   pfree(a);
01550   pfree(b);
01551 
01552   return(TRUE);
01553 }
01554 
01555 /* Contained */
01556 /* Box(A) Contained by Box(B) IFF Box(B) Contains Box(A) */
01557 bool geoobj_contained(GeoObj *a, GeoObj *b)
01558 {
01559 TRACE("geoobj_contained")
01560   if (geoobj_contains(b,a) == TRUE)
01561     return(TRUE);
01562   else
01563     return(FALSE);
01564 }
01565 
01566 /* Overlap */
01567 /* Box(A) Overlap Box(B) IFF (pt(a)LL < pt(B)UR) && (pt(b)LL < pt(a)UR) */
01568 bool geoobj_overlap(GeoObj *box_a, GeoObj *box_b)
01569 {
01570   int i;
01571   GeoObj *a;
01572   GeoObj *b;
01573   char *data;
01574   WKBGeometry *wkbgeoa;
01575   WKBGeometry *wkbgeob;
01576 
01577 TRACE("geoobj_overlap") 
01578 elog(NOTICE,geoobj_out(box_a));
01579 elog(NOTICE,geoobj_out(box_b));
01580  
01581   if ( (box_a==NULL) || (box_b==NULL) )
01582     return(FALSE);
01583 
01584 /*  a = swap_corners(box_a);
01585   b = swap_corners(box_b);
01586 */
01587 
01588   a=envelope(box_a);
01589   b=envelope(box_b);
01590 
01591   data=a->data;
01592   wkbgeoa=(WKBGeometry *)data;
01593 
01594   data=b->data;
01595   wkbgeob=(WKBGeometry *)data;
01596   
01597   /* compare within the dimensions of (b) */
01598   for ( i = 0; i < 3; i++ ) {
01599     if ( wkbgeoa->box.x[i] > wkbgeob->box.x[i + 3] )
01600       return(FALSE);
01601     if ( wkbgeoa->box.x[i + 3] < wkbgeob->box.x[i] )
01602       return(FALSE);
01603   }
01604 
01605   pfree(a);
01606   pfree(b);
01607 
01608   return(TRUE);
01609 }
01610 
01611 
01612 /* Distance */
01613 /* The distance is computed as a per axis sum of the squared distances
01614    between 1D projections of the boxes onto Cartesian axes. Assuming zero 
01615    distance between overlapping projections, this metric coincides with the 
01616    "common sense" geometric distance */
01617 float *geoobj_distance(GeoObj *a, GeoObj *b)
01618 {
01619   int i;
01620   double d, distance;
01621   float *result;
01622   char *data;
01623   WKBGeometry *wkbgeoa;
01624   WKBGeometry *wkbgeob;
01625 
01626 TRACE("geoobj_distance")
01627   data=a->data;
01628   wkbgeoa=(WKBGeometry *)data;
01629 
01630   data=b->data;
01631   wkbgeob=(WKBGeometry *)data;
01632   
01633   result = (float *) palloc(sizeof(float));
01634   
01635   distance = 0.0;
01636   /* compute within the dimensions of (b) */
01637   for ( i = 0; i < 3; i++ ) {
01638     d = distance_1D(wkbgeoa->box.x[i], wkbgeoa->box.x[i + 3], wkbgeob->box.x[i], wkbgeob->box.x[i + 3]);
01639     distance += d*d;
01640   }
01641 
01642   *result = (float)sqrt(distance);
01643 
01644   return(result);
01645 }
01646 
01647 static float distance_1D(float a1, float a2, float b1, float b2)
01648 {
01649 TRACE("distance_1D")
01650   /* interval (a) is entirely on the left of (b) */
01651   if( (a1 <= b1) && (a2 <= b1) && (a1 <= b2) && (a2 <= b2) ) {
01652     return ( min( b1, b2 ) - max( a1, a2 ) );
01653   }
01654 
01655   /* interval (a) is entirely on the right of (b) */
01656   if( (a1 > b1) && (a2 > b1) && (a1 > b2) && (a2 > b2) ) {
01657     return ( min( a1, a2 ) - max( b1, b2 ) );
01658   }
01659   
01660   /* the rest are all sorts of intersections */
01661   return(0.0);
01662 }
01663 
01664 /* normalize the box's co-ordinates by placing min(xLL,xUR) to LL
01665    and max(xLL,xUR) to UR 
01666 */
01667 static GeoObj *swap_corners( GeoObj *a )
01668 {
01669   int i, j;
01670   GeoObj * result;
01671   char *data;
01672   WKBGeometry *wkbgeoa;
01673   WKBGeometry *wkbgeoresult;
01674 
01675 TRACE("swap_corners")
01676   data=a->data;
01677   wkbgeoa=(WKBGeometry *)data;
01678     
01679   result = palloc(a->size);
01680   result->size = a->size;
01681   result->byte_order = BYTEORDER;
01682   result->wkbtype = FM_BOX;
01683   data=a->data;
01684   wkbgeoresult=(WKBGeometry *)data;
01685   wkbgeoresult->box.byteOrder = BYTEORDER;
01686   wkbgeoresult->box.wkbType = FM_BOX;
01687   
01688   for ( i = 0, j = 3; i < 3; i++, j++ ) {
01689     wkbgeoresult->box.x[i] = min(wkbgeoa->box.x[i],wkbgeoa->box.x[j]);
01690     wkbgeoresult->box.x[j] = max(wkbgeoa->box.x[i],wkbgeoa->box.x[j]);
01691   }
01692 
01693   return(result);
01694 }
01695 

Generated at Sat Jan 6 20:55:34 2001 for FMaps by doxygen1.2.1 written by Dimitri van Heesch, © 1997-2000