diff --git a/GeoIndex/GeoIndex.c b/GeoIndex/GeoIndex.c index a96feaf93e..5ac7a062d3 100644 --- a/GeoIndex/GeoIndex.c +++ b/GeoIndex/GeoIndex.c @@ -26,91 +26,62 @@ //////////////////////////////////////////////////////////////////////////////// /* GeoIndex.c - GeoIndex algorithms */ -/* Version 1.1 5.11.2011 R. A. Parker */ +/* Version 2.0 3.12.2011 R. A. Parker */ -#include "GeoIndex.h" +#ifdef GEO_TRIAGENS +#include "GeoIndex/GeoIndex.h" +#else +#include "GeoIndex2.h" +#endif /* Radius of the earth used for distances */ #define EARTHRADIUS 6371000.0 - /* geostring used to indicate end of free chain */ -#define GEOENDING -9000000000LL - - /* How big an "empty" index should be */ -#define GEOSTART 4 - /* How many results to allow for internally */ /* at the start of a "within distance" search */ /* (it grows if it needs to so is not critical */ -#define GEORESULTSTART 100 +#define GEORESULTSTART 100000 + +#define GEOSLOTSTART 50 +#define GEOPOTSTART 100 + +typedef struct +{ + GeoIndex * gi; + GeoCoordinate * gc; + double x; + double y; + double z; + GeoString gs; + GeoFix fixdist[GeoIndexFIXEDPOINTS]; + double snmd; + GeoFix distrej[GeoIndexFIXEDPOINTS]; +} GeoDetailedPoint; typedef struct { int pointsct; int allocpoints; - double maxdistance; int * slot; - double * distance; + double * snmd; } GeoResults; + +typedef struct +{ + GeoResults * gr; + GeoDetailedPoint * gd; + int stacksize; + int potid[46]; +} GeoStack; + + typedef struct { GeoIndex * gi; - GeoResults * gr; - int topofstack; - int startslotid[32]; - int endslotid[32]; - GeoString targetgs; - double targetlat; - double targetlong; -} GeoStack; - -/* Create a new GeoIndex for the user */ - -GeoIndex * GeoIndex_new(void) -{ - GeoIndex * gi; - int i; - - gi = malloc(sizeof(GeoIndex)); - if(gi==NULL) return gi; - gi->slotct = GEOSTART; - -/* try to allocate all the things we need */ - gi->sortslot = malloc(GEOSTART*sizeof(int)); - gi->geostrings = malloc(GEOSTART*sizeof(GeoString)); - gi->gc = malloc(GEOSTART*sizeof(GeoCoordinate)); - -/* if any of them fail, free the ones that succeeded */ -/* and then return the NULL pointer for our user */ - if ( ( gi->sortslot == NULL) || - ( gi->geostrings == NULL) || - ( gi->gc == NULL) ) - { - if ( gi->sortslot != NULL) free(gi->sortslot); - if ( gi->geostrings != NULL) free(gi->geostrings); - if ( gi->gc != NULL) free(gi->gc); - free(gi); - return NULL; - } - -/* initialize chain of empty slots */ - for(i=0;igeostrings[i]=-(i+1); - gi->geostrings[GEOSTART-1]=GEOENDING; - -/* set occslots to say there are no occupied slots */ -/* and therefore no used entries in the sortslots either */ - gi->occslots=0; - return gi; -} - -void GeoIndex_free(GeoIndex * gi) -{ - free(gi->gc); - free(gi->geostrings); - free(gi->sortslot); - free(gi); -} + int pathlength; + int path[46]; +} GeoPath; double GeoIndex_distance(GeoCoordinate * c1, GeoCoordinate * c2) { @@ -126,200 +97,284 @@ double GeoIndex_distance(GeoCoordinate * c1, GeoCoordinate * c2) return 2.0 * EARTHRADIUS * asin(mole/2.0); } -void GeoIndexFreeSlot(GeoIndex * gi, int slot) +void GeoIndexFreePot(GeoIndex * gi, int pot) { - gi->geostrings[slot]=gi->geostrings[0]; - gi->geostrings[0]= (GeoString) -slot; + gi->pots[pot].LorLeaf=gi->pots[0].LorLeaf; + gi->pots[0].LorLeaf = pot; } -int GeoIndexNewSlot(GeoIndex * gi) +int GeoIndexNewPot(GeoIndex * gi) { - long long x, y; - GeoString * gs; - GeoCoordinate * gc; - int * ss; - int newslotct,j; - if(gi->geostrings[0]==GEOENDING) + int newpotct,j; + long long x,y; + GeoPot * gp; + if( gi->pots[0].LorLeaf == 0) { /* do the growth calculation in long long to make sure it doesn't */ /* overflow when the size gets to be near 2^31 */ - x = gi->slotct; + x = gi->potct; y=100+GeoIndexGROW; x=x*y + 99; y=100; x=x/y; - if(x>2000000000L) return -2; - newslotct= (int) x; - gs = realloc(gi->geostrings,newslotct*sizeof(GeoString)); - gc = realloc(gi->gc,newslotct*sizeof(GeoCoordinate)); - ss = realloc(gi->sortslot,newslotct*sizeof(int)); - if(gs!=NULL) gi->geostrings=gs; - if(gc!=NULL) gi->gc=gc; - if(ss!=NULL) gi->sortslot=ss; - if(gs==NULL) return -2; - if(gc==NULL) return -2; - if(ss==NULL) return -2; - for(j=gi->slotct;jslotct=newslotct; + if(x>1000000000L) return -2; + newpotct= (int) x; + gp = realloc(gi->pots,newpotct*sizeof(GeoPot)); + if(gp!=NULL) gi->pots=gp; + else return -2; + for(j=gi->potct;jpotct=newpotct; } - j= (int) -(gi->geostrings[0]); - gi->geostrings[0]=gi->geostrings[j]; + j= gi->pots[0].LorLeaf; + gi->pots[0].LorLeaf=gi->pots[j].LorLeaf; return j; } +GeoIndex * GeoIndex_new(void) +{ + GeoIndex * gi; + int i,j; + double lat, lon, x, y, z; + + gi = malloc(sizeof(GeoIndex)); + if(gi==NULL) return gi; +/* try to allocate all the things we need */ + gi->pots = malloc(GEOPOTSTART*sizeof(GeoPot)); + gi->gc = malloc(GEOSLOTSTART*sizeof(GeoCoordinate)); + +/* if any of them fail, free the ones that succeeded */ +/* and then return the NULL pointer for our user */ + if ( ( gi->pots == NULL) || + ( gi->gc == NULL) ) + { + if ( gi->pots != NULL) free(gi->pots); + if ( gi->gc != NULL) free(gi->gc); + free(gi); + return NULL; + } + +/* initialize chain of empty slots */ + for(i=0;igc[i]).latitude=i+1; + else (gi->gc[i]).latitude=0; + } + +/* similarly set up free chain of empty pots */ + for(i=0;ipots[i].LorLeaf=i+1; + else gi->pots[i].LorLeaf=0; + } + + gi->potct = GEOPOTSTART; + gi->slotct = GEOSLOTSTART; + +/* set up the fixed points structure */ + +/* this code assumes GeoIndexFIXEDPOINTS is eight! */ + + for(i=0;ifixed.x)[i]=x; + (gi->fixed.y)[i]=y; + (gi->fixed.z)[i]=z; + } +/* set up the root pot */ + + j=GeoIndexNewPot(gi); + gi->pots[j].LorLeaf=0; /* leaf pot */ + gi->pots[j].RorPoints=0; /* with no points in it! */ + gi->pots[j].middle = 0ll; + gi->pots[j].start = 0ll; + gi->pots[j].end = 0x1FFFFFFFFFFFFFll; + gi->pots[j].level = 1; + for(i=0;ipots[j].maxdist[i]=0; + return gi; +} + +void GeoIndex_free(GeoIndex * gi) +{ + free(gi->gc); + free(gi->pots); + free(gi); +} + + /* 2^25 / 90 rounded down. Used to convert */ /* degrees of longitude and latitude into */ /* integers for use making a GeoString */ -#define STRINGPERDEGREE 372827.0222 +#define STRINGPERDEGREE 372827.01 +#define HILBERTMAX 67108863 +#define ARCSINFIX 41720.0 -GeoString GeoMkString(GeoCoordinate * c) +GeoString GeoMkHilbert(GeoCoordinate * c) { double x1,y1; - int x,y,i; - int xm,ym; GeoString z; - xm=0x4000000; - ym=0x2000000; - x1=c->longitude+180.0; - y1=c->latitude+90; + int x,y; + int i,nz,temp; + y1=c->latitude+90.0; + z=0; + x1=c->longitude; + if(c->longitude < 0.0) + { + x1=c->longitude+180.0; + z=1; + } x=(int) (x1*STRINGPERDEGREE); y=(int) (y1*STRINGPERDEGREE); - z=0; for(i=0;i<26;i++) { - z+= (x&xm) + (y&ym); - z=z<<1; - xm=xm>>1; - ym=ym>>1; - } - z+= (x&xm); - return z; -} - -void GeoStackSet (GeoIndex * gi, GeoStack * gk, - GeoCoordinate * c, GeoResults * gr) -{ - gk->gi = gi; - gk->gr = gr; - - if(gi->occslots==0) - { - gk->topofstack=-1; - } - else - { - gk->topofstack = 0; - gk->startslotid[0]=0; - gk->endslotid[0]=gi->occslots-1; - } - gk->targetgs = GeoMkString(c); - gk->targetlong = c->longitude; - gk->targetlat = c->latitude; -} - -/* Routine to look for the targetgs in the interval */ -/* that lies at the top of the stack - usually the */ -/* whole index, when the stack has just been set */ -/* The returned "int" is the highest sortslot such */ -/* that all higher ones are bigger GeoStrings */ -/* Note - this may be one less than the bottom of the */ -/* interval, and may indeed be -1 */ - -int GeoFind(GeoStack * gk) -{ - int sstop, ssbot, ssmid; - GeoIndex * gi; - if(gk->topofstack==-1) return -1; /* empty stack */ - ssbot=gk->startslotid[gk->topofstack]; - sstop=gk->endslotid[gk->topofstack]; - gi=gk->gi; - if(gk->targetgs < gi->geostrings[gi->sortslot[ssbot]]) - return ssbot-1; - if(gk->targetgs >= gi->geostrings[gi->sortslot[sstop]]) - return sstop; -/* from now on, ssbot is always too low, sstop always high enough */ - while(sstop>(ssbot+1)) - { - ssmid=(sstop+ssbot)/2; -/* if ssmid == sstop, it is ssbot that will be set to it */ - if(gk->targetgs < gi->geostrings[gi->sortslot[ssmid]]) - sstop=ssmid; - else - ssbot=ssmid; - } - return ssbot; -} - -int GeoIndex_insert(GeoIndex * gi, GeoCoordinate * c) -{ - GeoStack gk; - int topslot,putslot,slot; - if(c->longitude < -180.0) return -3; - if(c->longitude > 180.0) return -3; - if(c->latitude < -90.0) return -3; - if(c->latitude > 90.0) return -3; - GeoStackSet(gi, &gk, c, (GeoResults *)NULL); - topslot = GeoFind(&gk); - putslot = topslot+1; /* we plan to put it here */ - -/* there may be any number of equal entries here, so we must */ -/* test them to ensure that we do not do a duplicate insertion */ - while(topslot > -1) - { - slot=gi->sortslot[topslot]; - if(gi->geostrings[slot] < gk.targetgs) break; - if( (c->latitude == ((gi->gc)[slot]).latitude) && - (c->longitude == ((gi->gc)[slot]).longitude) && - (c->data == ((gi->gc)[slot]).data) ) return -1; - topslot--; - } -/* All OK. we're going to try to put it in. */ - slot=GeoIndexNewSlot(gi); - if(slot==-2) return -2; /* no room :( */ - gi->geostrings[slot]=gk.targetgs; - (gi->gc)[slot].latitude = c->latitude; - (gi->gc)[slot].longitude=c->longitude; - (gi->gc)[slot].data=c->data; - -/* shuffle up the sortslots to make room */ - memmove(gi->sortslot+putslot+1, - gi->sortslot+putslot, - (gi->occslots-putslot)*sizeof(int)); - -/* put it in */ - gi->sortslot[putslot]=slot; - gi->occslots++; - return 0; -} - -int GeoIndex_remove(GeoIndex * gi, GeoCoordinate * c) -{ - int sortslot,slot; - GeoStack gk; - GeoStackSet(gi, &gk, c, (GeoResults *) NULL); - sortslot = GeoFind(&gk); - -/* there may be any number of equal entries here, so we must */ -/* test them all to find the one we want */ - while(sortslot > -1) - { - slot=gi->sortslot[sortslot]; - if(gi->geostrings[slot]>gk.targetgs) return -1; - if( (c->latitude == (gi->gc)[slot].latitude) && - (c->longitude == (gi->gc)[slot].longitude) && - (c->data == (gi->gc)[slot].data) ) + z<<=2; + nz=((y>>24)&2)+(x>>25); + x = (x<<1)&(HILBERTMAX); + y = (y<<1)&(HILBERTMAX); + if(nz==0) { - GeoIndexFreeSlot(gi,slot); - memmove(gi->sortslot+sortslot, - gi->sortslot+sortslot+1, - (gi->occslots-sortslot-1)*sizeof(int)); - - gi->occslots--; - return 0; + temp=x; + x=y; + y=temp; + } + if(nz==1) + { + temp=HILBERTMAX-x; + x=HILBERTMAX-y; + y=temp; + z+=3; + } + if(nz==2) + { + z+=1; + } + if(nz==3) + { + z+=2; } - sortslot--; } - return -1; + return z+1ll; + +} + +void GeoMkDetail(GeoIndex * gi, GeoDetailedPoint * gd, GeoCoordinate * c) +{ +/* entire routine takes about 0.94 microseconds */ + double x1,y1,z1,snmd; + int i; + gd->gi=gi; + gd->gc = c; +/* The GeoString computation takes about 0.17 microseconds */ + gd->gs=GeoMkHilbert(c); +/* This part takes about 0.32 microseconds */ + gd->z=sin(c->latitude*M_PI/180.0); + gd->x=cos(c->latitude*M_PI/180.0)*cos(c->longitude*M_PI/180.0); + gd->y=cos(c->latitude*M_PI/180.0)*sin(c->longitude*M_PI/180.0); +/* and this bit takes about 0.45 microseconds */ + for(i=0;ifixed.x)[i]; + y1=(gi->fixed.y)[i]; + z1=(gi->fixed.z)[i]; + snmd=(x1-gd->x)*(x1-gd->x)+(y1-gd->y)*(y1-gd->y)+ + (z1-gd->z)*(z1-gd->z); + (gd->fixdist)[i] = asin(sqrt(snmd)/2.0)*ARCSINFIX; + } +} + +double GeoMetersToSNMD(double meters) +{ + double angle,hnmd; + angle=0.5*meters/EARTHRADIUS; + hnmd=sin(angle); /* half normalized mole distance */ + if(angle>=M_PI/2.0) return 4.0; + else return hnmd*hnmd*4.0; +} + +void GeoSetDistance(GeoDetailedPoint * gd, double snmd) +{ + GeoFix gf; + int i; + gd->snmd = snmd; + gf = asin(sqrt(snmd)/2.0)*ARCSINFIX; + gf++; + for(i=0;ifixdist)[i]<=gf) (gd->distrej)[i]=0; + else (gd->distrej)[i]=(gd->fixdist)[i]-gf; + } +} + +void GeoStackSet (GeoStack * gk, GeoDetailedPoint * gd, GeoResults * gr) +{ + int pot; + GeoIndex * gi; + GeoPot * gp; + gi=gd->gi; + gk->gr = gr; + gk->gd = gd; + gk->stacksize = 0; + pot=1; + while(1) + { + gp=gi->pots+pot; + if(gp->LorLeaf==0) break; + if(gp->middle>gd->gs) + { + gk->potid[gk->stacksize]=gp->RorPoints; + pot=gp->LorLeaf; + } + else + { + gk->potid[gk->stacksize]=gp->LorLeaf; + pot=gp->RorPoints; + } + gk->stacksize++; + } + gk->potid[gk->stacksize]=pot; } GeoResults * GeoResultsCons(int alloc) @@ -342,11 +397,70 @@ GeoResults * GeoResultsCons(int alloc) gres->pointsct = 0; gres->allocpoints = alloc; gres->slot = sa; - gres->distance = dd; -/* no need to initialize maxdistance */ + gres->snmd = dd; +/* no need to initialize maxsnmd */ return gres; } +void GeoResultsStartCount(GeoResults * gr) +{ + int i; + for(i=0;iallocpoints;i++) + { + gr->slot[i]=0; + gr->snmd[i]=10.0; + } +} + +void GeoResultsInsertPoint(GeoResults * gr, int slot, double snmd) +{ + int i,j1,j2,temp; + if(snmd>=gr->snmd[0]) return; + if(gr->slot[0]==0) gr->pointsct++; + i=0; /* i is now considered empty */ + while(1) + { + j1=2*i+1; + j2=2*i+2; + if(j1allocpoints) + { + if(j2allocpoints) + { + if(gr->snmd[j1]>gr->snmd[j2]) + { + temp=j1; + j1=j2; + j2=temp; + } + /* so now j2 is >= j1 */ + if(gr->snmd[j2]<=snmd) + { + gr->snmd[i]=snmd; + gr->slot[i]=slot; + return; + } + gr->snmd[i]=gr->snmd[j2]; + gr->slot[i]=gr->slot[j2]; + i=j2; + continue; + } + if(gr->snmd[j1]<=snmd) + { + gr->snmd[i]=snmd; + gr->slot[i]=slot; + return; + } + gr->snmd[i]=gr->snmd[j1]; + gr->slot[i]=gr->slot[j1]; + i=j1; + continue; + } + gr->snmd[i]=snmd; + gr->slot[i]=slot; + return; + } +} + int GeoResultsGrow(GeoResults * gr) { int newsiz; @@ -357,15 +471,15 @@ int GeoResultsGrow(GeoResults * gr) newsiz=gr->pointsct + (gr->pointsct/2) + 1; if(newsiz > 1000000000) return -1; sa=realloc(gr->slot, newsiz*sizeof(int)); - dd=realloc(gr->distance, newsiz*sizeof(double)); + dd=realloc(gr->snmd, newsiz*sizeof(double)); if( (sa==NULL) || (dd==NULL) ) { if(sa!=NULL) gr->slot = sa; - if(dd!=NULL) gr->distance = dd; + if(dd!=NULL) gr->snmd = dd; return -1; } gr->slot = sa; - gr->distance = dd; + gr->snmd = dd; gr->allocpoints = newsiz; return 0; } @@ -374,7 +488,8 @@ GeoCoordinates * GeoAnswers (GeoIndex * gi, GeoResults * gr) { GeoCoordinates * ans; GeoCoordinate * gc; - int i; + int i,j,slot; + double mole; ans = malloc(sizeof(GeoCoordinates));; gc = malloc(gr->pointsct * sizeof(GeoCoordinate)); if( (ans==NULL) || (gc==NULL) ) @@ -382,194 +497,104 @@ GeoCoordinates * GeoAnswers (GeoIndex * gi, GeoResults * gr) if(ans!=NULL) free(ans); if(gc!=NULL) free(gc); free(gr->slot); - free(gr->distance); + free(gr->snmd); free(gr); return NULL; } ans->length = gr->pointsct; ans->coordinates = gc; - ans->distances = gr->distance; - for(i=0;ipointsct;i++) + j=0; + for(i=0;iallocpoints;i++) { - ans->coordinates[i].latitude = - (gi->gc)[gr->slot[i]].latitude; - ans->coordinates[i].longitude = - (gi->gc)[gr->slot[i]].longitude; - ans->coordinates[i].data = - (gi->gc)[gr->slot[i]].data; + if(j>=gr->pointsct) break; + slot=gr->slot[i]; + if(slot==0) continue; + ans->coordinates[j].latitude = + (gi->gc)[slot].latitude; + ans->coordinates[j].longitude = + (gi->gc)[slot].longitude; + ans->coordinates[j].data = + (gi->gc)[slot].data; + mole=sqrt(gr->snmd[i]); + if(mole > 2.0) mole = 2.0; /* make sure arcsin succeeds! */ + gr->snmd[j]= 2.0 * EARTHRADIUS * asin(mole/2.0); + j++; } + ans->distances = gr->snmd; free(gr->slot); free(gr); return ans; } -/* "Reverse-Engineer" a GeoString - convert it into */ -/* the given GeoCoordinate */ - -void GeoReverse(GeoString gs, GeoCoordinate * gc) +int GeoPotJunk(GeoDetailedPoint * gd, int pot) { - int lon,lat; - int i,j,k; - double d; - j=0; - lon=0; - lat=0; - -/* un-interleave the bits of the GeoString */ - - for(i=0;i<27;i++) - { - k=(gs>>j)&3; - lon+=(k&1)<>=1; - lat+=(k&1)<longitude=(d/STRINGPERDEGREE)-180.0; - d=lat; - gc->latitude=(d/STRINGPERDEGREE)-90.0; + int i; + GeoPot * gp; + gp=(gd->gi)->pots + pot; + for(i=0;imaxdist[i]distrej[i]) return 1; + return 0; } -double GeoMinDistance(GeoStack * gk, GeoString k1, GeoString k2) +double GeoSNMD(GeoDetailedPoint * gd, GeoCoordinate * c) { - GeoString xor,msb,mst; - int bits; - GeoCoordinate d1, d2; - - double md,md1; - double latdiff; - double lon1, lon2; - double x,y,phi; - GeoCoordinate c1, c2; -/* if the point is in the interval, fast answer clearly zero */ - if( (k1<=gk->targetgs) && (k2>=gk->targetgs) ) return 0.0; - -/* compute mst - (top) mask for bits where k1 and k2 agree */ -/* msb - (bottom) mask for lower bits where they may differ */ - - xor = k1 ^ k2; - bits=0; - while(xor!=0) - { - bits++; - xor>>=1; - } - msb=(0xFFFFFFFFFFFFFFLL)>>(56-bits); - mst=(0xFFFFFFFFFFFFFFLL)^msb; -/* compute d1 as coordinates of lower corner */ - GeoReverse(k1&mst,&d1); -/* and d2 as coordinates of upper corner */ - GeoReverse((k1&mst)+msb,&d2); -/* d1 negative d2 positive values */ - if( (d1.longitude<=gk->targetlong) && (d2.longitude>=gk->targetlong) ) - { - if(d1.latitude<=gk->targetlat) - { - if(d2.latitude>=gk->targetlat) - return 0.0; /* in the box */ - latdiff=d2.latitude-gk->targetlat; - } - else latdiff=d1.latitude-gk->targetlat; - md = EARTHRADIUS*M_PI*fabs(latdiff)/180.0 - 3.0; - if(md<0.0) md=0.0; - return md; - } - lon1=fabs(d1.longitude-gk->targetlong); - lon2=fabs(d2.longitude-gk->targetlong); - if(lon1>180.0) lon1=360.0-lon1; - if(lon2>180.0) lon2=360.0-lon2; - if(lon2targetlat; - x=cos(gk->targetlat*M_PI/180.0)*cos(lon1*M_PI/180.0); - y=sin(gk->targetlat*M_PI/180.0); - phi=atan2(y,x)*180.0/M_PI; - if( phi>=d1.latitude ) - { - if(phi<=d2.latitude) c2.latitude=phi; - else c2.latitude=d2.latitude; - } - else c2.latitude=d1.latitude; - md=GeoIndex_distance(&c1,&c2) - 2.0; - if(md<0.0) md=0.0; - return md; - } - c1.longitude=0.0; - c2.longitude=lon1; - c2.latitude=gk->targetlat; - c1.latitude=d1.latitude; - md=GeoIndex_distance(&c1,&c2) - 2.0; - c1.latitude=d2.latitude; - md1=GeoIndex_distance(&c1,&c2) - 2.0; - if(md1latitude*M_PI/180.0); + x=cos(c->latitude*M_PI/180.0)*cos(c->longitude*M_PI/180.0); + y=cos(c->latitude*M_PI/180.0)*sin(c->longitude*M_PI/180.0); + return (x-gd->x)*(x-gd->x) + (y-gd->y)*(y-gd->y) + + (z-gd->z)*(z-gd->z); } GeoCoordinates * GeoIndex_PointsWithinRadius(GeoIndex * gi, GeoCoordinate * c, double d) { - GeoStack gk; - GeoCoordinates * answer; GeoResults * gres; - double dist; - int i,slot,r; - int sstop, ssbot, ssmid; - GeoString gstop, gsbot; + GeoCoordinates * answer; + GeoDetailedPoint gd; + GeoStack gk; + GeoPot * gp; + int r,pot,slot,i; + double snmd,maxsnmd; gres=GeoResultsCons(100); if(gres==NULL) return NULL; - GeoStackSet(gi,&gk,c,gres); - while (gk.topofstack > -1) + GeoMkDetail(gi,&gd,c); + GeoStackSet(&gk,&gd,gres); + maxsnmd=GeoMetersToSNMD(d); + GeoSetDistance(&gd,maxsnmd); + gk.stacksize++; + while(gk.stacksize>=1) { - ssbot=gk.startslotid[gk.topofstack]; - gsbot=gi->geostrings[gi->sortslot[ssbot]]; - sstop=gk.endslotid[gk.topofstack]; - gstop=gi->geostrings[gi->sortslot[sstop]]; - dist=GeoMinDistance(&gk,gsbot,gstop); - if(dist>d) - { - gk.topofstack--; + gk.stacksize--; + pot=gk.potid[gk.stacksize]; + if(GeoPotJunk(&gd,pot)) continue; - } - if( (sstop-ssbot) < GeoIndexSMALLINTERVAL ) + gp=gi->pots+pot; + if(gp->LorLeaf==0) { - for(i=ssbot;i<=sstop;i++) + for(i=0;iRorPoints;i++) { - slot=gi->sortslot[i]; - dist=GeoIndex_distance(c,gi->gc + slot); - if(dist<=d) + slot=gp->points[i]; + snmd=GeoSNMD(&gd,gi->gc+slot); + if(snmd>maxsnmd) continue; + r = GeoResultsGrow(gres); + if(r==-1) { - r=GeoResultsGrow(gres); - if(r==-1) - { - free(gres->distance); - free(gres->slot); - free(gres); - return NULL; - } - gres->slot[gres->pointsct]=slot; - gres->distance[gres->pointsct]=dist; - gres->pointsct++; + free(gres->snmd); + free(gres->slot); + free(gres); + return NULL; } + gres->slot[gres->pointsct]=slot; + gres->snmd[gres->pointsct]=snmd; + gres->pointsct++; } - gk.topofstack--; - continue; } - ssmid=(sstop+ssbot)/2; - gk.startslotid[gk.topofstack]=ssbot; - gk.endslotid[gk.topofstack]=ssmid; - - gk.topofstack++; - gk.startslotid[gk.topofstack]=ssmid+1; - gk.endslotid[gk.topofstack]=sstop; - + else + { + gk.potid[gk.stacksize++]=gp->LorLeaf; + gk.potid[gk.stacksize++]=gp->RorPoints; + } } answer=GeoAnswers(gi,gres); return answer; /* note - this may be NULL */ @@ -578,114 +603,461 @@ GeoCoordinates * GeoIndex_PointsWithinRadius(GeoIndex * gi, GeoCoordinates * GeoIndex_NearestCountPoints(GeoIndex * gi, GeoCoordinate * c, int count) { - GeoResults * gres; + GeoResults * gr; + GeoDetailedPoint gd; GeoCoordinates * answer; GeoStack gk; - GeoString gsbot, gstop; - int ssbot, sstop, ssmid, slot, i, k; - double maxdist, dist, newmaxdist; - gres=GeoResultsCons(count); - if(gres==NULL) return NULL; - GeoStackSet(gi,&gk,c,gres); -/* Phase 1 - get initial stack with nearer intervals maybe at top */ - while (gk.topofstack > -1) /*allowing for an empty index */ - { - ssbot=gk.startslotid[gk.topofstack]; - sstop=gk.endslotid[gk.topofstack]; - if( (sstop-ssbot) < GeoIndexSMALLINTERVAL ) break; - ssmid=(sstop+ssbot)/2; - slot=gi->sortslot[ssmid]; - if(gk.targetgs > gi->geostrings[slot]) - { - gk.startslotid[gk.topofstack]=ssbot; - gk.endslotid[gk.topofstack]=ssmid; - gk.topofstack++; - gk.startslotid[gk.topofstack]=ssmid+1; - gk.endslotid[gk.topofstack]=sstop; - } - else - { - gk.startslotid[gk.topofstack]=ssmid+1; - gk.endslotid[gk.topofstack]=sstop; - gk.topofstack++; - gk.startslotid[gk.topofstack]=ssbot; - gk.endslotid[gk.topofstack]=ssmid; - } - } + GeoPot * gp; + int pot,slot,i,left; + double snmd; -/* Phase 2 - populate the results list with first points found */ - maxdist=0.0; /* maximum distance = low-values */ - while (gk.topofstack > -1) /* allowing for early emptying */ + gr=GeoResultsCons(count); + if(gr==NULL) return NULL; + GeoMkDetail(gi,&gd,c); + GeoStackSet(&gk,&gd,gr); + GeoResultsStartCount(gr); + left=count; + + while(gk.stacksize>=0) { - ssbot=gk.startslotid[gk.topofstack]; - sstop=gk.endslotid[gk.topofstack]; - slot=gi->sortslot[ssbot]; - dist=GeoIndex_distance(c,gi->gc + slot); - if(dist>maxdist) maxdist=dist; - gres->slot[gres->pointsct]=slot; - gres->distance[gres->pointsct]=dist; - gres->pointsct++; - ssbot++; - if(ssbot>sstop) + pot=gk.potid[gk.stacksize--]; + gp=gi->pots+pot; + if(left<=0) { - gk.topofstack--; + GeoSetDistance(&gd,gr->snmd[0]); + if(GeoPotJunk(&gd,pot)) continue; + } + if(gp->LorLeaf==0) + { + for(i=0;iRorPoints;i++) + { + slot=gp->points[i]; + snmd=GeoSNMD(&gd,gi->gc+slot); + GeoResultsInsertPoint(gr,slot,snmd); + left--; + if(left<-1) left=-1; + } } else { - gk.startslotid[gk.topofstack]=ssbot; - } - if(gres->pointsct>=count) break; /* normal termination */ - } -/* Phase 3 - put in the rest of the points as needed */ -/* maxdist is the largest value in the list, so we must beat that */ - while (gk.topofstack > -1) /* done when stack empty */ - { - ssbot=gk.startslotid[gk.topofstack]; - gsbot=gi->geostrings[gi->sortslot[ssbot]]; - sstop=gk.endslotid[gk.topofstack]; - gstop=gi->geostrings[gi->sortslot[sstop]]; - dist=GeoMinDistance(&gk,gsbot,gstop); - if(dist>=maxdist) - { - gk.topofstack--; /* throw away the whole interval */ - continue; - } - if( (sstop-ssbot) < GeoIndexSMALLINTERVAL ) - { - for(i=ssbot;i<=sstop;i++) + if(gd.gs>gp->middle) { - slot=gi->sortslot[i]; - dist=GeoIndex_distance(c,gi->gc + slot); - if(dist<=maxdist) + gk.potid[++gk.stacksize]=gp->LorLeaf; + gk.potid[++gk.stacksize]=gp->RorPoints; + } + else + { + gk.potid[++gk.stacksize]=gp->RorPoints; + gk.potid[++gk.stacksize]=gp->LorLeaf; + } + } + } + answer=GeoAnswers(gi,gr); + return answer; /* note - this may be NULL */ +} + +void GeoIndexFreeSlot(GeoIndex * gi, int slot) +{ + gi->gc[slot].latitude=gi->gc[0].latitude; + gi->gc[0].latitude = slot; +} + +int GeoIndexNewSlot(GeoIndex * gi) +{ + int newslotct,j; + long long x,y; + GeoCoordinate * gc; + if(gi->gc[0].latitude == 0.0) + { +/* do the growth calculation in long long to make sure it doesn't */ +/* overflow when the size gets to be near 2^31 */ + x = gi->slotct; + y=100+GeoIndexGROW; + x=x*y + 99; + y=100; + x=x/y; + if(x>2000000000L) return -2; + newslotct= (int) x; + gc = realloc(gi->gc,newslotct*sizeof(GeoCoordinate)); + if(gc!=NULL) gi->gc=gc; + else return -2; + for(j=gi->slotct;jslotct=newslotct; + } + j= (int) (gi->gc[0].latitude); + gi->gc[0].latitude=gi->gc[j].latitude; + return j; +} + +int GeoFind(GeoPath * gt, GeoDetailedPoint * gd) +{ + int pot,pot1; + int i; + int slot; + GeoIndex * gi; + GeoCoordinate * gc; + GeoPot * gp; + gi = gd->gi; + gt->gi = gi; + pot=1; + gt->pathlength=0; + while(1) + { + gp=gi->pots+pot; + gt->path[gt->pathlength]=pot; + gt->pathlength++; + if(gp->LorLeaf == 0) break; + if(gp->middle > gd->gs) + pot=gp->LorLeaf; + else + pot=gp->RorPoints; + } +/* so we have a pot such that top is bigger but bottom isn't */ + while(1) /* so look for an exact match */ + { + for(i=0;iRorPoints;i++) + { + slot=gp->points[i]; + gc=gi->gc+slot; + if( ( (gd->gc)->latitude ==gc->latitude ) && + ( (gd->gc)->longitude==gc->longitude) && + ( (gd->gc)->data ==gc->data ) ) + { + gt->path[gt->pathlength]=i; + return 1; + } + } + if(gp->startgs) break; +/* need to find the predecessor of this pot */ +/* this is expected to be a rare event, so */ +/* no time is wasted to simplify this! */ + while(1) + { + gt->pathlength--; + pot1=gt->path[gt->pathlength-1]; + gp=gi->pots+pot1; + if(pot==gp->RorPoints) break; /* cannot go off the front */ + pot=pot1; + } + gp=gi->pots+pot1; + pot=gp->LorLeaf; +/* now we have a pot whose iterated right child we want */ + while(1) + { + gp=gi->pots+pot; + gt->path[gt->pathlength]=pot; + gt->pathlength++; + if(gp->LorLeaf == 0) break; + pot=gp->RorPoints; + } + } + return 2; +} + +void GeoPopulateMaxdist(GeoIndex * gi, GeoPot * gp, GeoString * gsa) +{ + int i,j; + GeoDetailedPoint gd; + gsa[0]=0x1FFFFFFFFFFFFFll; + gsa[1]=0ll; + for(j=0;jmaxdist[j]=0; + for(i=0;iRorPoints;i++) + { + GeoMkDetail(gi,&gd,gi->gc+gp->points[i]); + for(j=0;jgp->maxdist[j]) + gp->maxdist[j] = gd.fixdist[j]; + if(gd.gsgsa[1]) gsa[1]=gd.gs; + } + gp->level=1; +} + +int GeoGetPot(GeoPath * gt, int height) +{ + return gt->path[gt->pathlength-height]; +} + +void GeoAdjust(GeoIndex * gi, int potx) /* the kids are alright */ +{ + int poty,potz; /* x = (yz) */ + int i; + GeoPot * gpx; + GeoPot * gpy; + GeoPot * gpz; + gpx=gi->pots + potx; + poty=gpx->LorLeaf; + gpy=gi->pots + poty; + potz=gpx->RorPoints; + gpz=gi->pots + potz; + gpx->start=gpy->start; + gpx->end =gpz->end; + gpx->middle=gpz->start; + gpx->level=gpy->level; + if( (gpz->level) > gpx->level) + gpx->level = gpz->level; + gpx->level++; + for(i=0;imaxdist[i]=gpy->maxdist[i]; + if( gpx->maxdist[i]maxdist[i]) + gpx->maxdist[i]=gpz->maxdist[i]; + + } +} + +void RotateLeft(GeoIndex * gi, int pote) +{ + int pota,potb,potc,potd; + GeoPot * gpd; + GeoPot * gpe; + gpe=gi->pots + pote; + potd=gpe->RorPoints; + gpd=gi->pots + potd; + pota=gpe->LorLeaf; + potb=gpd->LorLeaf; + potc=gpd->RorPoints; + gpd->LorLeaf=pota; + gpd->RorPoints=potb; + GeoAdjust(gi,potd); + gpe->LorLeaf=potd; + gpe->RorPoints=potc; + GeoAdjust(gi,pote); +} + +void RotateRight(GeoIndex * gi, int pote) +{ + int pota,potb,potc,potd; + GeoPot * gpd; + GeoPot * gpe; + gpe=gi->pots + pote; + potd=gpe->LorLeaf; + gpd=gi->pots + potd; + pota=gpd->LorLeaf; + potb=gpd->RorPoints; + potc=gpe->RorPoints; + gpd->LorLeaf=potb; + gpd->RorPoints=potc; + GeoAdjust(gi,potd); + gpe->LorLeaf=pota; + gpe->RorPoints=potd; + GeoAdjust(gi,pote); +} + +int GeoIndex_insert(GeoIndex * gi, GeoCoordinate * c) +{ + int i,j,js,slot,pot,pot1,pot2; + int potx,pota,poty,potz; + int lvx,lv1,lva,lvy,lvz; + int height,rebalance; + GeoDetailedPoint gd; + GeoPath gt; + GeoPot * gp; + GeoPot * gp1; + GeoPot * gp2; + GeoPot * gpx; + GeoPot * gpy; + GeoPot * gpz; + GeoPot * gpa; + GeoString gsa[2]; + GeoString mings, gs; + rebalance=0; + if(c->longitude < -180.0) return -3; + if(c->longitude > 180.0) return -3; + if(c->latitude < -90.0) return -3; + if(c->latitude > 90.0) return -3; + GeoMkDetail(gi,&gd,c); + i = GeoFind(>,&gd); + if(i==1) return -1; + pot=gt.path[gt.pathlength-1]; + gp=gi->pots + pot; +/* new point, so we try to put it in */ + slot=GeoIndexNewSlot(gi); + if(slot==-2) return -2; /* no room :( */ + gi->gc[slot].latitude =c->latitude; + gi->gc[slot].longitude=c->longitude; + gi->gc[slot].data =c->data; +/* check first if we are going to need two new pots, and */ +/* if we are, go get them now before we get too tangled */ + if(gp->RorPoints == GeoIndexPOTSIZE) + { + rebalance=1; + pot1=GeoIndexNewPot(gi); + pot2=GeoIndexNewPot(gi); + gp=gi->pots + pot; /* may have re-alloced! */ + if( (pot1==-2) || (pot2==-2) ) + { + GeoIndexFreeSlot(gi,slot); + if(pot1!=-2) GeoIndexFreePot(gi,pot1); + if(pot2!=-2) GeoIndexFreePot(gi,pot2); + return -2; + } + gp1=gi->pots + pot1; + gp2=gi->pots + pot2; +/* pot is old one, pot1 and pot2 are the new ones */ + gp1->LorLeaf=0; /* leaf pot */ + gp1->RorPoints=0; /* no points in it yet */ +/* first move the points from pot to pot2 */ + gp2->LorLeaf=0; /* leaf pot */ + gp2->RorPoints=gp->RorPoints; + for(i=0;iRorPoints;i++) + gp2->points[i]=gp->points[i]; +/* move the first half of the points from pot2 to pot1 */ + for(i=0;iRorPoints;j++) + { + gs=GeoMkHilbert(gi->gc + gp2->points[j]); + if(gsdistance[k] == maxdist) - { - gres->slot[k]=slot; - gres->distance[k]=dist; - dist=1000000000.0; /* put in just once */ - } - if(gres->distance[k]>newmaxdist) - newmaxdist=gres->distance[k]; - } - maxdist=newmaxdist; + mings=gs; + js=j; } } - gk.topofstack--; - continue; + gp1->points[gp1->RorPoints]=gp2->points[js]; + gp2->points[js]=gp2->points[gp2->RorPoints-1]; + gp2->RorPoints--; + gp1->RorPoints++; + } + GeoPopulateMaxdist(gi,gp2,gsa); + mings=gsa[0]; + GeoPopulateMaxdist(gi,gp1,gsa); + mings=(mings+gsa[1])/2ll; + gp1->start=gp->start; + gp1->end=mings; + gp2->start=mings; + gp2->end=gp->end; + gp->LorLeaf=pot1; + gp->RorPoints=pot2; + GeoAdjust(gi,pot); + gt.pathlength++; + if(gd.gspoints[gp->RorPoints]=slot; + gp->RorPoints++; +/* now propagate the maxdistances */ + for(i=0;i=0) + { + if(gd.fixdist[i] > gi->pots[gt.path[j]].maxdist[i]) + gi->pots[gt.path[j]].maxdist[i] = gd.fixdist[i]; + else break; + j--; + } + } +/* just need to balance the tree */ + if(rebalance==0) return 0; + height=2; + while(1) + { + potx=GeoGetPot(>,height); + gpx=gi->pots + potx; + lvx=gpx->level; + if(potx==1) break; + /* root pot ? */ + pot1=GeoGetPot(>,height+1); /* pot1=parent(x) */ + gp1=gi->pots + pot1; + lv1=gp1->level; + if(lv1>lvx) break; + if(gp1->LorLeaf==potx) /* gpx is the left child? */ + { + pota=gp1->RorPoints; /* 1 = (xa) */ + gpa=gi->pots+pota; + lva=gpa->level; + if( (lva+1) == lv1) /* so it is legal to up lev(1) */ + { + gp1->level++; + height++; + continue; + } + poty=gpx->RorPoints; + gpy=gi->pots + poty; + lvy=gpy->level; + potz=gpx->LorLeaf; + gpz=gi->pots + potz; + lvz=gpz->level; + if(lvy<=lvz) + { + RotateRight(gi,pot1); + height++; + continue; + } + RotateLeft(gi,potx); + RotateRight(gi,pot1); + } + else /* gpx is the right child */ + { + pota=gp1->LorLeaf; /* 1 = (ax) */ + gpa=gi->pots+pota; + lva=gpa->level; + if( (lva+1) == lv1) /* so it is legal to up lev(1) */ + { + gp1->level++; + height++; + continue; + } + poty=gpx->LorLeaf; + gpy=gi->pots + poty; + lvy=gpy->level; + potz=gpx->RorPoints; + gpz=gi->pots + potz; + lvz=gpz->level; + if(lvy<=lvz) + { + RotateLeft(gi,pot1); + height++; + continue; + } + RotateRight(gi,potx); + RotateLeft(gi,pot1); + } + } + return 0; +} + +int GeoIndex_remove(GeoIndex * gi, GeoCoordinate * c) +{ + GeoDetailedPoint gd; + GeoPot * gp; + GeoPath gt; + int i,pot,potix,slot,pathix; + GeoMkDetail(gi,&gd,c); + i = GeoFind(>,&gd); + if(i!=1) return -1; + pot=gt.path[gt.pathlength-1]; + gp=gi->pots + pot; + potix = gt.path[gt.pathlength]; + slot=gp->points[potix]; + GeoIndexFreeSlot(gi,slot); + gp->points[potix]=gp->points[gp->RorPoints-1]; + gp->RorPoints--; + if( (2*gp->RorPoints)0) + { + pathix--; + pot=gt.path[pathix]; + GeoAdjust(gi,pot); + } + return 0; } void GeoIndex_CoordinatesFree(GeoCoordinates * clist) @@ -695,65 +1067,61 @@ void GeoIndex_CoordinatesFree(GeoCoordinates * clist) free(clist); } -/* the following code can all be deleted if required */ +int GeoIndex_hint(GeoIndex * gi, int hint) +{ + return 0; +} #ifdef DEBUG +void RecursivePotDump (GeoIndex * gi, FILE * f, int pot) +{ + int i; + GeoPot * gp; + GeoCoordinate * gc; + gp=gi->pots + pot; + fprintf(f,"GP. pot %d level %d Kids %d %d\n", + pot, gp->level, gp->LorLeaf,gp->RorPoints); + fprintf(f,"strings %llx %llx %llx\n",gp->start,gp->middle,gp->end); + fprintf(f,"maxdists "); + for(i=0;imaxdist[i]); + fprintf(f,"\n"); + if(gp->LorLeaf==0) + { + fprintf(f,"Leaf pot containing %d points . . .\n",gp->RorPoints); + for(i=0;iRorPoints;i++) + { + fprintf(f,"Child %d Point %d ",i,gp->points[i]); + gc=gi->gc + gp->points[i]; + fprintf(f,"Lat. %9.4f, Long. %9.4f", + gc->latitude, gc->longitude); +#if DEBUG==2 + fprintf(f," %s",(char *) gc->data); +#endif + fprintf(f,"\n"); + } + } + else + { + fprintf(f,"\nPot %d - Left Child of pot %d\n", + gp->LorLeaf,pot); + RecursivePotDump (gi,f,gp->LorLeaf); + fprintf(f,"\nPot %d - Right Child of pot %d\n", + gp->RorPoints,pot); + RecursivePotDump (gi,f,gp->RorPoints); + } +} + void GeoIndex_INDEXDUMP (GeoIndex * gi, FILE * f) { - int i,j,k,dig; - GeoString gs; - fprintf(f,"Maximum Slots currently %d\n",gi->slotct); - fprintf(f,"Number of entries in the index = %d,",gi->occslots); - fprintf(f,"\n"); - fprintf(f,"Latitude Longitude GeoString "); -#if DEBUG == 2 - fprintf(f," Data"); -#endif - fprintf(f,"\n"); - for(i=0;ioccslots;i++) - { - j=gi->sortslot[i]; - fprintf(f,"%7.2f %7.2f ", - gi->gc[j].latitude,gi->gc[j].longitude); - gs=gi->geostrings[j]; - for(k=0;k<18;k++) - { - dig=(gs>>58)&15; - fprintf(f,"%x",dig); - gs<<=4; - } -#if DEBUG == 2 - fprintf(f," %s",(char *)gi->gc[j].data); -#endif - fprintf(f,"\n"); - } + fprintf(f,"Dump of entire index. %d pots and %d slots allocated\n", + gi->potct, gi->slotct); + RecursivePotDump(gi,f,1); } int GeoIndex_INDEXVALID(GeoIndex * gi) { - int i,j; - GeoString gs,freeslot; - for(i=0;ioccslots;i++) - { - j=gi->sortslot[i]; - gs = GeoMkString(&(gi->gc[j])); - if(gi->geostrings[j]!=gs) - return -1; - } - gs=0; - for(i=0;ioccslots;i++) - { - j=gi->sortslot[i]; - if(gi->geostrings[j]geostrings[j]; - } - freeslot=0; - for(i=gi->occslots+1;islotct;i++) - freeslot=-gi->geostrings[freeslot]; - if(gi->geostrings[freeslot]!=GEOENDING) - return -3; return 0; } diff --git a/GeoIndex/GeoIndex.h b/GeoIndex/GeoIndex.h index 9379c78556..f02e7d0c92 100644 --- a/GeoIndex/GeoIndex.h +++ b/GeoIndex/GeoIndex.h @@ -26,7 +26,7 @@ //////////////////////////////////////////////////////////////////////////////// /* GeoIndex.h - header file for GeoIndex algorithms */ -/* Version 1.0 24.9.2011 R. A. Parker */ +/* Version 2.0 3.12.2011 R. A. Parker */ #ifdef GEO_TRIAGENS #include @@ -43,10 +43,14 @@ typedef long long GeoString; /* percentage growth of slot or slotslot tables */ -#define GeoIndexGROW 50 +#define GeoIndexGROW 10 -/* Largest interval that is processed simply */ -#define GeoIndexSMALLINTERVAL 5 +/* maximum number of points in a pot >=2 */ +#define GeoIndexPOTSIZE 6 + +/* number of fixed points for max-dist calculations */ +#define GeoIndexFIXEDPOINTS 8 +typedef unsigned short GeoFix; /* If this #define is there, then the INDEXDUMP and */ /* INDEXVALID functions are also available. These */ @@ -55,13 +59,13 @@ typedef long long GeoString; /* assumed to be a character string, if DEBUG is */ /* set to 2. */ -/* #define DEBUG 1 */ +#define DEBUG 1 typedef struct { double latitude; double longitude; - void const * data; + void * data; } GeoCoordinate; typedef struct @@ -73,11 +77,29 @@ typedef struct typedef struct { - int slotct; - int occslots; /* number of occupied slots */ - int * sortslot; - GeoString * geostrings; /* These two indexed in */ - GeoCoordinate * gc; /* parallel by a "slot" */ + double x[GeoIndexFIXEDPOINTS]; + double y[GeoIndexFIXEDPOINTS]; + double z[GeoIndexFIXEDPOINTS]; +} GeoIndexFixed; + +typedef struct +{ + int LorLeaf; + int RorPoints; + GeoString middle; + GeoFix maxdist[GeoIndexFIXEDPOINTS]; + GeoString start; + GeoString end; + int level; + int points[GeoIndexPOTSIZE]; +} GeoPot; +typedef struct +{ + GeoIndexFixed fixed; /* fixed point data */ + int potct; /* pots allocated */ + int slotct; /* slots allocated */ + GeoPot * pots; /* the pots themselves */ + GeoCoordinate * gc; /* the slots themselves */ } GeoIndex; GeoIndex * GeoIndex_new(void); @@ -85,6 +107,7 @@ void GeoIndex_free(GeoIndex * gi); double GeoIndex_distance(GeoCoordinate * c1, GeoCoordinate * c2); int GeoIndex_insert(GeoIndex * gi, GeoCoordinate * c); int GeoIndex_remove(GeoIndex * gi, GeoCoordinate * c); +int GeoIndex_hint(GeoIndex * gi, int hint); GeoCoordinates * GeoIndex_PointsWithinRadius(GeoIndex * gi, GeoCoordinate * c, double d); GeoCoordinates * GeoIndex_NearestCountPoints(GeoIndex * gi, diff --git a/UnitTests/georeg.c b/UnitTests/georeg.c index 28778d88ac..cc8942a7e9 100644 --- a/UnitTests/georeg.c +++ b/UnitTests/georeg.c @@ -26,12 +26,12 @@ //////////////////////////////////////////////////////////////////////////////// /* regression testing program for GeoIndex module */ -/* R.A.P. 1.1 5.11.2011 */ +/* R.A.P. 2.0 4.12.2011 */ #ifdef GEO_TRIAGENS #include "GeoIndex/GeoIndex.h" #else -#include "GeoIndex.h" +#include "GeoIndex2.h" #endif int errors; @@ -39,6 +39,24 @@ int errors; char ix[1000]; /* working array for (void *) data */ char iy[30]; +int np[20]={14676,23724,24457,24785,24809, + 25909,25051,27052,25192,28126, + 25191,29291,25079,30338,24852, + 31320,24565,31936,24133,31893}; +int hs[20]={16278595,6622245,83009659,97313687,94638085, + 86998272,72133969,99595673,76554853,116384443, + 458789,67013205,44378533,97387502,97331554, + 76392514,43104144,48695421,87361440,1556675} ; +int np1[10]={6761,13521,20282,27042,33803, + 40563,47324,54084,60845,67605}; +int hs1[10]={79121168,71376992,120173779,54504134,89075457, + 50229904,41454125,104395668,14196041,106196305}; +int hs2[10]={21216619,99404510,92771863,40046216,25926851, + 6056147,93877377,82650316,14776130,41666384}; +int np4[4]={2838,5116,5180,9869}; +int hs4[4]={33972992,9770664,11661062,28398735}; +int hs5[4]={79685116,67516870,19274248,35037618}; + void icheck(int e, int a, int b) { if(a==b) return; @@ -118,6 +136,56 @@ void gcmass(int e, GeoCoordinates * gc, int ct, int hash) GeoIndex_CoordinatesFree(gc); } +void coonum(GeoCoordinate * gc, int num) +{ + double lat,lon; + int i,j; + j=num; + gc->latitude=-42.23994323; + gc->longitude=-53.40029372; + gc->data=&ix[num%997]; + for(i=1;i<30;i++) + { + lat=0.0; + lon=0.0; + if( (j&1)==1 ) + { + if(i==1) { lat=16.1543722; lon=12.1992384; } + if(i==2) { lat=6.14347227; lon=5.19923843; } + if(i==3) { lat=2.15237222; lon=3.19923843; } + if(i==4) { lat=1.51233226; lon=1.63122143; } + if(i==5) { lat=1.14347229; lon=1.69932121; } + if(i==6) { lat=0.93443431; lon=0.80023323; } + if(i==7) { lat=0.63443442; lon=0.79932032; } + if(i==8) { lat=0.70049323; lon=0.80076994; } + if(i==9) { lat=0.89223236; lon=0.78805238; } + if(i==10) { lat=0.85595656; lon=0.72332312; } + if(i==11) { lat=0.18823232; lon=0.21129576; } + if(i==12) { lat=0.12294854; lon=0.15543207; } + if(i==13) { lat=0.30984343; lon=0.18223745; } + if(i==14) { lat=0.19923412; lon=0.27345381; } + if(i==15) { lat=0.18223534; lon=0.15332087; } + if(i==16) { lat=0.09344343; lon=0.08002332; } + if(i==17) { lat=0.06344344; lon=0.07993203; } + if(i==18) { lat=0.07004932; lon=0.08007699; } + if(i==19) { lat=0.08922323; lon=0.07880523; } + if(i==20) { lat=0.08559565; lon=0.07233231; } + if(i==21) { lat=0.01882323; lon=0.02112957; } + if(i==22) { lat=0.01229485; lon=0.01554320; } + if(i==23) { lat=0.03098434; lon=0.01822374; } + if(i==24) { lat=0.01992341; lon=0.02734538; } + if(i==25) { lat=0.01822353; lon=0.01533208; } + if(i==26) { lat=0.00934434; lon=0.00800233; } + if(i==27) { lat=0.00634434; lon=0.00799320; } + if(i==28) { lat=0.00700493; lon=0.00800769; } + if(i==29) { lat=0.00480493; lon=0.00170769; } + } + gc->latitude+=lat; + gc->longitude+=lon; + j>>=1; + } +} + int main(int argc, char ** argv) { GeoIndex * gi; @@ -170,7 +238,7 @@ int main(int argc, char ** argv) /* mainly for debugging rather than regression */ gi=GeoIndex_new(); - + i=GeoIndex_hint(gi,10); /* set it to "robust" mode */ for(i=0;i<50;i++) { gcp.latitude = 90.0; @@ -900,44 +968,27 @@ int main(int argc, char ** argv) if(mode==2) { - gcp1.latitude = 44.999999; - gcp1.longitude= 45.0000001; - printf("start of 50000 x points within radius (43 found)\n"); - for(i=0;i<50000;i++) + printf("start of 5000000 x points within radius 127 Km (3 found)\n"); + for(i=0;i<10000;i++) { - list1 = GeoIndex_PointsWithinRadius(gi,&gcp1,430000.0); - GeoIndex_CoordinatesFree(list1); + for(j=1;j<=500;j++) + { + coonum(&gcp1,j); + list1 = GeoIndex_PointsWithinRadius(gi,&gcp1,127000.0); + GeoIndex_CoordinatesFree(list1); + } } printf("End of timing test\n"); } if(mode==3) { - printf("start of 10000 x points by count (2)\n"); - for(i=0;i<4;i++) + printf("start of 5000000 x points by count (2)\n"); + for(i=0;i<10000;i++) { - if(i==0) - { - gcp1.latitude = 44.999999; - gcp1.longitude= 45.0000001; - } - if(i==1) - { - gcp1.latitude = -44.999999; - gcp1.longitude= 45.0000001; - } - if(i==2) - { - gcp1.latitude = 44.999999; - gcp1.longitude= -45.0000001; - } - if(i==3) - { - gcp1.latitude = -44.999999; - gcp1.longitude= -45.0000001; - } - for(j=0;j<2500;j++) + for(j=1;j<=500;j++) { + coonum(&gcp1,j); list1 = GeoIndex_NearestCountPoints(gi,&gcp1,2); GeoIndex_CoordinatesFree(list1); } @@ -949,6 +1000,134 @@ int main(int argc, char ** argv) GeoIndex_free(gi); +/* */ +/* 400 - 499 */ +/* ========= */ +/* */ + +/* another batch of massive tests - again used as */ +/* comparisons - this time between 1.1 and 2.0 */ + + gi=GeoIndex_new(); + la=41.23456789; + lo=39.87654321; + for(i=0;i<20;i++) + { + for(j=1;j<30000;j++) + { + gcp.latitude = la; + gcp.longitude= lo; + gcp.data = &ix[(7*i+j)%1000]; + r = GeoIndex_insert(gi,&gcp); + icheck(400,0,r); + la+=19.5396157761; + if(la>90.0) la-=180.0; + lo+=17.2329155421; + if(lo>180) lo-=360.0; + } + list1 = GeoIndex_PointsWithinRadius(gi,&gcp,9800000.0); + for(j=0;jlength;j++) /* delete before freeing list1! */ + { + r=GeoIndex_remove(gi,list1->coordinates+j); + } + gcmass(400+5*i,list1, np[i], hs[i]); + } + + GeoIndex_free(gi); + +/* */ +/* 500 - 599 */ +/* ========= */ +/* */ + +/* This set of tests aims to cluster the points in */ +/* a sligtly more realistic way. */ + + gi=GeoIndex_new(); + for(i=1;i<135212;i++) + { + coonum(&gcp,i); + r = GeoIndex_insert(gi,&gcp); + icheck(501,0,r); + if( (i%2)==0) + { + coonum(&gcp,i/2); + r = GeoIndex_remove(gi,&gcp); + icheck(502,0,r); + } + if( (i%13521)==0) + { + list1 = GeoIndex_PointsWithinRadius(gi,&gcp,9800000.0); + gcmass(505+4*(i/13521),list1, + np1[i/13521 - 1], hs1[i/13521 - 1]); + list1 = GeoIndex_NearestCountPoints(gi,&gcp,i/1703); + gcmass(507+4*(i/13521),list1,i/1703,hs2[i/13521 - 1]); + } + } + + GeoIndex_free(gi); +/* */ +/* 600 - 610 */ +/* ========= */ +/* */ + +/* Test a very large database - too big for V1.2! */ + if(mode==4) + { + gi=GeoIndex_new(); + for(i=1;i<21123212;i++) + { + coonum(&gcp,i); + r = GeoIndex_insert(gi,&gcp); + icheck(601,0,r); + if( (i%2)==0) + { + coonum(&gcp,i/2); + r = GeoIndex_remove(gi,&gcp); + icheck(602,0,r); + } + if( (i%4541672)==0) + { + list1 = GeoIndex_PointsWithinRadius(gi,&gcp,9800.0); + gcmass(603+4*(i/4541672),list1, + np4[i/4541672 - 1], hs4[i/4541672 - 1]); + list1 = GeoIndex_NearestCountPoints(gi,&gcp,i/1832703); + gcmass(605+4*(i/4541672),list1,i/1832703,hs5[i/4541672 - 1]); + } + } + GeoIndex_free(gi); + } +/* */ +/* 900 - 999 */ +/* ========= */ +/* */ + +/* tests that ensure old bugs have not reappeared */ + +/* forgot to allow for distance greater than pi.radius */ +/* this went into sin(theta) to give small value, so it */ +/* found that many points were not within 30000 Km */ + + gi=GeoIndex_new(); + la=41.23456789; + lo=39.87654321; + for(j=1;j<50;j++) + { + gcp.latitude = la; + gcp.longitude= lo; + gcp.data = &ix[j]; + r = GeoIndex_insert(gi,&gcp); + icheck(900,0,r); + la+=19.5396157761; + if(la>90.0) la-=180.0; + lo+=17.2329155421; + if(lo>180) lo-=360.0; + } + list1 = GeoIndex_PointsWithinRadius(gi,&gcp,30000000.0); + gcmass(901,list1, 49, 94065911); + + GeoIndex_free(gi); + if(errors==0) { printf("Georeg Regression Program succeeded\n"); return 0;