Hawkin
/
arangodb
mirror of https://gitee.com/bigwinds/arangodb
1
0
Fork 0
Code Issues Projects Releases Wiki Activity
arangodb/arangod/GeoIndex/GeoIndex.cpp

2294 lines
80 KiB
C++
Raw Blame History

////////////////////////////////////////////////////////////////////////////////
/// DISCLAIMER
///
/// Copyright 2014-2016 ArangoDB GmbH, Cologne, Germany
/// Copyright 2004-2014 triAGENS GmbH, Cologne, Germany
///
/// Licensed under the Apache License, Version 2.0 (the "License");
/// you may not use this file except in compliance with the License.
/// You may obtain a copy of the License at
///
/// http://www.apache.org/licenses/LICENSE-2.0
///
/// Unless required by applicable law or agreed to in writing, software
/// distributed under the License is distributed on an "AS IS" BASIS,
/// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
/// See the License for the specific language governing permissions and
/// limitations under the License.
///
/// Copyright holder is ArangoDB GmbH, Cologne, Germany
///
/// @author R. A. Parker
////////////////////////////////////////////////////////////////////////////////
/* GeoIndex.c - GeoIndex algorithms */
/* Version 2.1 8.1.2012 R. A. Parker */
#define _USE_MATH_DEFINES
#include <math.h>
#include "GeoIndex.h"
/* Radius of the earth used for distances */
#define EARTHRADIAN 6371000.0
#define GEOSLOTSTART 50
#define GEOPOTSTART 100
#if GeoIndexFIXEDSET == 2
#define GeoIndexFIXEDPOINTS 2
#endif
#if GeoIndexFIXEDSET == 3
#define GeoIndexFIXEDPOINTS 3
#endif
#if GeoIndexFIXEDSET == 4
#define GeoIndexFIXEDPOINTS 4
#endif
#if GeoIndexFIXEDSET == 5
#define GeoIndexFIXEDPOINTS 5
#endif
#if GeoIndexFIXEDSET == 6
#define GeoIndexFIXEDPOINTS 6
#endif
#if GeoIndexFIXEDSET == 8
#define GeoIndexFIXEDPOINTS 8
#endif
#ifndef GeoIndexFIXEDPOINTS
#define GeoIndexFIXEDPOINTS 1
#endif
/* =================================================== */
/* GeoIndexFixed structure. */
/* Only occurs once, and that is in the GeoIx struct */
/* holds the x,y and z coordinates (between -1 and +1) */
/* of the fixed points used for pot rejection purposes */
/* They are computed at GeoIndex_new time and not */
/* changed after that */
/* =================================================== */
typedef struct {
double x[GeoIndexFIXEDPOINTS];
double y[GeoIndexFIXEDPOINTS];
double z[GeoIndexFIXEDPOINTS];
} GeoIndexFixed;
/* =================================================== */
/* GeoPot structure */
/* These only occur in the main index itself, and the */
/* GeoIx structure has an array of them. The data */
/* items are arranged so that the access during a */
/* search is approximately sequential, which should be */
/* a little faster on most machines. */
/* The first two data items are used for several */
/* different purposes. LorLeaf is zero for a leaf pot */
/* and the left child for a non-leaf pot. RorPoints */
/* is the right child for a non-leaf pot, and the */
/* number of points in the pot for a leaf pot */
/* The three GeoString values give the bounds (weak) */
/* for the Hilbert values in this pot. middle is not */
/* used for a leaf pot. */
/* maxdist is the maximum, over all points descendent */
/* from this pot, of the distances to the fixed points */
/* level is the AVL-level. It is 1 for a leaf pot, */
/* and always at least 1 more and at most 2 more than */
/* each of its children, and exactly 1 more than at */
/* least one of its children, - the AVL spec. */
/* "points" lists the slotid of the points. This is */
/* only used for a leaf pot. */
/* =================================================== */
typedef struct {
int LorLeaf;
int RorPoints;
GeoString middle;
GeoFix maxdist[GeoIndexFIXEDPOINTS];
GeoString start;
GeoString end;
int level;
int points[GeoIndexPOTSIZE];
} GeoPot;
/* =================================================== */
/* GeoIx structure */
/* This is the REAL GeoIndex structure - the one in */
/* the GeoIndex.h file is just a sham (it says it is */
/* a char!) to keep the structure private so that the */
/* GeoIndex.h is short and contains only data of */
/* interest to the user. */
/* The GeoIx structure basically consists of two */
/* arrays - the slots (the points) and the pots (the */
/* Balanced (AVL) search tree for finding near points) */
/* The Fixed-point data is held here also, giving the */
/* x, y and z coordinates of the fixed points, this */
/* data being the fastest to use */
/* potct and slotct are used when the index needs to */
/* grow (because it has run out of slots or pots) */
/* There is no provision at present for the index to */
/* get smaller when the majority of points are deleted */
/* =================================================== */
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 */
size_t _memoryUsed; /* the amount of memory currently used */
} GeoIx;
/* =================================================== */
/* GeoDetailedPoint structure */
/* The routine GeoMkDetail is given a point - really */
/* just a latitude and longitude, and computes all the */
/* values in this GeoDetailedPoint structure. */
/* This is intended to include everything that will be */
/* needed about the point, and is called both for the */
/* searches (count and distance) and the updates */
/* (insert and remove). It is only ever useful */
/* locally - it is created, populated, used and */
/* forgotten all within a single user's call */
/* the GeoIx is noted there to simplify some calls */
/* The GeoCoordinate (a pointer to the user's one) */
/* is included. The x, y and z coordinates (between */
/* 1 and -1) are computed, as is the GeoString - the */
/* Hilbert curve value used to decide where in the */
/* index a point belongs. The fixdist array is the */
/* distance to the fixed points. */
/* The other two entries (snmd and distrej) are not */
/* computed by GeoMkDetail, but are put put in place */
/* later, for the searches only, by GeoSetDistance. */
/* They basically hold the radius of the circle around */
/* the target point outside which indexed points will */
/* be too far to be of interest. This is set once and */
/* for all in the case of a search-by-distance, but */
/* for a search-by-count the interesting distance */
/* decreases as further points are found. */
/* Anyway, snmd hold the radius in SNMD form (squared */
/* normalized mole distance) being the distance in */
/* three-dimensional space between two points passing */
/* through the earth (as a mole digs!) - this being */
/* the fastest to compute on the fly, and is used for */
/* looking at individual points to decide whether to */
/* include them. The distrej array, on the other hand */
/* is the array of distances to the fixed points, and */
/* is used to reject pots (leaf or non-leaf). */
/* The routine GeoPotJunk is used to test this, */
/* by comparing the distances in the pot the this array*/
/* =================================================== */
typedef struct {
GeoIx* gix;
GeoCoordinate* gc;
double x;
double y;
double z;
GeoString gs;
GeoFix fixdist[GeoIndexFIXEDPOINTS];
double snmd;
GeoFix distrej[GeoIndexFIXEDPOINTS];
} GeoDetailedPoint;
/* =================================================== */
/* GeoResults structure */
/* During the searches, this structure is used to */
/* accumulate the points that will be returned */
/* In the case of a search-by-distance, the results are*/
/* simply a list, which is grown by about 50% if the */
/* initial allocation (100) is inadequte. In the case */
/* of a search-by-count, the exact number needed is */
/* known from the start, but the structure is not just */
/* a simple list in this case. Instead it is organized*/
/* as a "priority queue" to enable large values of the */
/* <count> parameter to be rapidly processed. In the */
/* case of count, each value is kept to be larger that */
/* both of its "children" - at 2n+1 and 2n+2. Hence */
/* the largest distance is always at position 0 and can*/
/* be readily found, but if it is to be replaced, there*/
/* is some procession (no more than log(count) work) */
/* to do to find the correct place to insert the new */
/* one in the priority queue. This work is done in the*/
/* GeoResultsInsertPoint routine (not used by distance)*/
/* =================================================== */
typedef struct {
int pointsct;
int allocpoints;
int* slot;
double* snmd;
} GeoResults;
/* =================================================== */
/* GeoStack structure */
/* During searches of both kinds, at any time there is */
/* this "stack" (first-in-last-out) of pots still to be*/
/* processed. At the start of a search of either type,*/
/* this structure is populated (by GeoStackSet) by */
/* starting at the root pot, and selecting a child that*/
/* could contain the target point. The other pot is */
/* put on the stack and processing continues. The */
/* stack is then processed by taking a pot off, */
/* discarding it if the maximum distance to a fixed */
/* point is too low, and otherwise putting both the */
/* children onto the stack (since it is faster to do */
/* this than suffer the cache miss to determine whether*/
/* either or both of the children can be rejected) */
/* =================================================== */
typedef struct {
GeoResults* gr;
GeoDetailedPoint* gd;
int stacksize;
int potid[50];
} GeoStack;
/* =================================================== */
/* GeoPath structure */
/* Similar in many ways to the GeoStack, above, this */
/* structure is used during insertion and deletion. */
/* Notice that the pots of the index to not contain */
/* pointers to their parent, since this is not needed */
/* during a search. During insertion and removal, */
/* however, it is necessary to move upwards to */
/* propogate the maximum distances and to balance the */
/* tree. Hence the GeoFind procedure, called at the */
/* beginning of insertion and deletion, populates this */
/* structure so that the full path from the root node */
/* to the current pot being considered is known, and */
/* its parent found when needed. */
/* =================================================== */
typedef struct {
GeoIx* gix;
int pathlength;
int path[50];
} GeoPath;
/* =================================================== */
/* GeoIndex_Distance routine */
/* This is the user-facing routine to compute the */
/* distance in meters between any two points, given */
/* by latitude and longitude in a pair of GeoCoordinate*/
/* structures. It operates by first converting the */
/* two points into x, y and z coordinates in 3-space, */
/* then computing the distance between them (again in */
/* three space) using Pythagoras, computing the angle */
/* subtended at the earth's centre, between the two */
/* points, and finally muliply this angle (in radians) */
/* by the earth's radius to convert it into meters. */
/* =================================================== */
double GeoIndex_distance(GeoCoordinate* c1, GeoCoordinate* c2) {
/* math.h under MacOS defines y1 and j1 as global variable */
double xx1, yy1, z1, x2, y2, z2, mole;
z1 = sin(c1->latitude * M_PI / 180.0);
xx1 = cos(c1->latitude * M_PI / 180.0) * cos(c1->longitude * M_PI / 180.0);
yy1 = cos(c1->latitude * M_PI / 180.0) * sin(c1->longitude * M_PI / 180.0);
z2 = sin(c2->latitude * M_PI / 180.0);
x2 = cos(c2->latitude * M_PI / 180.0) * cos(c2->longitude * M_PI / 180.0);
y2 = cos(c2->latitude * M_PI / 180.0) * sin(c2->longitude * M_PI / 180.0);
mole = sqrt((xx1 - x2) * (xx1 - x2) + (yy1 - y2) * (yy1 - y2) +
(z1 - z2) * (z1 - z2));
if (mole > 2.0) mole = 2.0; /* make sure arcsin succeeds! */
return 2.0 * EARTHRADIAN * asin(mole / 2.0);
}
/* =================================================== */
/* GeoIndexFreePot */
/* takes the supplied pot, and puts it back onto the */
/* free list. */
/* =================================================== */
void GeoIndexFreePot(GeoIx* gix, int pot) {
gix->pots[pot].LorLeaf = gix->pots[0].LorLeaf;
gix->pots[0].LorLeaf = pot;
}
/* =================================================== */
/* GeoIndexNewPot */
/* During insertion, it may happen that a leaf pot */
/* becomes full. In this case this routine is called */
/* (always twice, as it happens) to allocate a new */
/* leaf pot, and a new pot to become the parent of both*/
/* the old and the new leaf pots. Usually this will */
/* be a simple matter of taking a pot off the free */
/* list, but occasionally the free list will be empty, */
/* in which case the pot array must be realloced. */
/* NOTICE that in this case, the pots may have moved, */
/* so it is critically important ot ensure that any */
/* pointers to pots are re-computed after this routine */
/* has been called! The GeoIndex_insert routine is */
/* therefore careful to get the new pots (if any are */
/* needed) before it gets too far into things. */
/* =================================================== */
int GeoIndexNewPot(GeoIx* gix) {
int j;
GeoPot* gp;
if (gix->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 */
long long x = gix->potct;
long long y = 100 + GeoIndexGROW;
x = x * y + 99;
y = 100;
x = x / y;
if (x > 1000000000L) return -2;
int newpotct = (int)x;
gp = static_cast<GeoPot*>(TRI_Reallocate(TRI_UNKNOWN_MEM_ZONE, gix->pots,
newpotct * sizeof(GeoPot)));
if (gp == NULL) {
return -2;
}
gix->pots = gp;
// update memory usage
gix->_memoryUsed -= gix->potct * sizeof(GeoPot);
gix->_memoryUsed += newpotct * sizeof(GeoPot);
for (j = gix->potct; j < newpotct; j++) {
GeoIndexFreePot(gix, j);
}
gix->potct = newpotct;
}
j = gix->pots[0].LorLeaf;
gix->pots[0].LorLeaf = gix->pots[j].LorLeaf;
return j;
}
/* =================================================== */
/* GeoIndex_new routine */
/* User-facing routine to create a whole new GeoIndex. */
/* Much of the bulk of the code in this routine is */
/* populating the fixed points, depending on which */
/* set of fixed points are in used. */
/* The first job is to allocate the initial arrays for */
/* holding the points, and the pots that index them. */
/* If this fails, no harm is done and the NULL pointer */
/* is returned. Otherwise all the point and pots are */
/* put onto their respective free lists. */
/* The fixed point structure is then set up. */
/* Finally the root pot (pot 1) is set up to be a leaf */
/* pot containing no points, but with the start and end*/
/* GeoString values (points on the Hilbert Curve) set */
/* to be "low values" and "high values" respectively, */
/* being slightly outside the range of possible */
/* GeoString values of real (latitude, longitude) */
/* points */
/* =================================================== */
GeoIdx* GeoIndex_new(void) {
GeoIx* gix;
int i, j;
gix = static_cast<GeoIx*>(
TRI_Allocate(TRI_UNKNOWN_MEM_ZONE, sizeof(GeoIx), false));
if (gix == NULL) {
return (GeoIdx*)gix;
}
/* try to allocate all the things we need */
gix->pots = static_cast<GeoPot*>(
TRI_Allocate(TRI_UNKNOWN_MEM_ZONE, GEOPOTSTART * sizeof(GeoPot), false));
gix->gc = static_cast<GeoCoordinate*>(TRI_Allocate(
TRI_UNKNOWN_MEM_ZONE, GEOSLOTSTART * sizeof(GeoCoordinate), false));
/* if any of them fail, free the ones that succeeded */
/* and then return the NULL pointer for our user */
if ((gix->pots == NULL) || (gix->gc == NULL)) {
if (gix->pots != NULL) {
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gix->pots);
}
if (gix->gc != NULL) {
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gix->gc);
}
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gix);
return NULL;
}
// set initial memory usage
gix->_memoryUsed =
GEOPOTSTART * sizeof(GeoPot) + GEOSLOTSTART * sizeof(GeoCoordinate);
/* initialize chain of empty slots */
for (i = 0; i < GEOSLOTSTART; i++) {
if (i < GEOSLOTSTART - 1)
(gix->gc[i]).latitude = i + 1;
else
(gix->gc[i]).latitude = 0;
}
/* similarly set up free chain of empty pots */
for (i = 0; i < GEOPOTSTART; i++) {
if (i < GEOPOTSTART - 1)
gix->pots[i].LorLeaf = i + 1;
else
gix->pots[i].LorLeaf = 0;
}
gix->potct = GEOPOTSTART;
gix->slotct = GEOSLOTSTART;
/* set up the fixed points structure */
for (i = 0; i < GeoIndexFIXEDPOINTS; i++) {
double lat = 90.0;
double lon = 0.0;
#if GeoIndexFIXEDSET == 2
if (i == 1) {
lat = -90.0;
lon = 0.0;
}
#endif
#if GeoIndexFIXEDSET == 3
if (i == 1) {
lat = -30.0;
lon = 0.0;
}
if (i == 2) {
lat = -30;
lon = 180.0;
}
#endif
#if GeoIndexFIXEDSET == 4
if (i == 1) {
lat = -19.471220634490691369246;
lon = 180.0;
}
if (i == 2) {
lat = -19.471220634490691369246;
lon = -60.0;
}
if (i == 3) {
lat = -19.471220634490691369246;
lon = 60.0;
}
#endif
#if GeoIndexFIXEDSET == 5
if (i == 1) {
lat = -90.0;
lon = 0.0;
}
if (i == 2) {
lat = 0.0;
lon = 0.0;
}
if (i == 3) {
lat = 0.0;
lon = 120.0;
}
if (i == 4) {
lat = 0.0;
lon = -120.0;
}
#endif
#if GeoIndexFIXEDSET == 6
if (i == 1) {
lat = -90.0;
lon = 0.0;
}
if (i == 2) {
lat = 0.0;
lon = 0.0;
}
if (i == 3) {
lat = 0.0;
lon = 180.0;
}
if (i == 4) {
lat = 0.0;
lon = 90.0;
}
if (i == 5) {
lat = 0.0;
lon = -90.0;
}
#endif
#if GeoIndexFIXEDSET == 8
if (i == 1) {
lat = -90.0;
lon = 0.0;
}
if (i == 2) {
lat = 19.471220634490691369246;
lon = 0.0;
}
if (i == 3) {
lat = -19.471220634490691369246;
lon = 180.0;
}
if (i == 4) {
lat = 19.471220634490691369246;
lon = 120.0;
}
if (i == 5) {
lat = -19.471220634490691369246;
lon = -60.0;
}
if (i == 6) {
lat = 19.471220634490691369246;
lon = -120.0;
}
if (i == 7) {
lat = -19.471220634490691369246;
lon = 60.0;
}
#endif
double z = sin(lat * M_PI / 180.0);
double x = cos(lat * M_PI / 180.0) * cos(lon * M_PI / 180.0);
double y = cos(lat * M_PI / 180.0) * sin(lon * M_PI / 180.0);
(gix->fixed.x)[i] = x;
(gix->fixed.y)[i] = y;
(gix->fixed.z)[i] = z;
}
/* set up the root pot */
j = GeoIndexNewPot(gix);
gix->pots[j].LorLeaf = 0; /* leaf pot */
gix->pots[j].RorPoints = 0; /* with no points in it! */
gix->pots[j].middle = 0ll;
gix->pots[j].start = 0ll;
gix->pots[j].end = 0x1FFFFFFFFFFFFFll;
gix->pots[j].level = 1;
for (i = 0; i < GeoIndexFIXEDPOINTS; i++) gix->pots[j].maxdist[i] = 0;
return (GeoIdx*)gix;
}
/* =================================================== */
/* GeoIndex_free routine */
/* Destroys the GeoIndex, and frees all the memory that*/
/* this GeoIndex system allocated. Note that any */
/* objects that may have been pointed to by the user's */
/* data pointers are (of course) not freed by this call*/
/* =================================================== */
void GeoIndex_free(GeoIdx* gi) {
GeoIx* gix;
if (gi == NULL) {
return;
}
gix = (GeoIx*)gi;
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gix->gc);
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gix->pots);
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gix);
}
/* =================================================== */
/* GeoMkHilbert routine */
/* Points in this system are indexed by the "GeoString */
/* value, which is the distance to the point along the */
/* Hilbert Curve. This space-filling curve is best */
/* understood in a square, where the curve joins the */
/* bottom left to the bottom right. It consists of */
/* four copies of the Hilbert curve, one in each of the*/
/* four squares, going via the points half-way up the */
/* left side, the middle of the (large) square and half*/
/* way up the right side. Notice that the first and */
/* last of these are flipped on a diagonal, whereas the*/
/* middle two, going along the top half, are in the */
/* original orientation, but at half the size. This */
/* description matches the code below, except that the */
/* two hemispheres are imagined to be squares where the*/
/* poles are the top line and the bottom line of the */
/* square. */
/* =================================================== */
/* 2^25 / 90 rounded down. Used to convert */
/* degrees of longitude and latitude into */
/* integers for use making a GeoString */
#define STRINGPERDEGREE 372827.01
/* 2^26 - 1 = 0x3ffffff */
#define HILBERTMAX 67108863
GeoString GeoMkHilbert(GeoCoordinate* c) {
/* math.h under MacOS defines y1 and j1 as global variable */
double xx1, yy1;
GeoString z;
int x, y;
int i, temp;
yy1 = c->latitude + 90.0;
z = 0;
xx1 = c->longitude;
if (c->longitude < 0.0) {
xx1 = c->longitude + 180.0;
z = 1;
}
x = (int)(xx1 * STRINGPERDEGREE);
y = (int)(yy1 * STRINGPERDEGREE);
for (i = 0; i < 26; i++) {
z <<= 2;
int nz = ((y >> 24) & 2) + (x >> 25);
x = (x << 1) & (HILBERTMAX);
y = (y << 1) & (HILBERTMAX);
if (nz == 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;
}
}
return z + 1ll;
}
/* =================================================== */
/* GeoMkDetail routine */
/* At the beginning of both searches, and also at the */
/* start of an insert or remove, this routine is called*/
/* to compute all the detail that can usefully be found*/
/* once and for all. The timings below were on done on*/
/* a 2011 ordinary desktop pentium */
/* 0.94 microseconds is - very approximately - 20% of */
/* the execution time of searches and/or updates, so */
/* is an obvious target for future speedups should they*/
/* be required (possibly by using less-accurate trig. */
/* it consists of three essentially separate tasks */
/* 1. Find the GeoString (Hilbert) value. */
/* 2. compute the x, y and z coordinates */
/* 3. find the distances to the fixed points */
/* all of these are needed for all of the operations */
/* =================================================== */
#if GEOFIXLEN == 2
#define ARCSINFIX 41720.0
/* resolution about 300 meters */
#endif
#if GEOFIXLEN == 4
#define ARCSINFIX 1520000000.0
/* resolution about 3 cm */
#endif
void GeoMkDetail(GeoIx* gix, GeoDetailedPoint* gd, GeoCoordinate* c) {
/* entire routine takes about 0.94 microseconds */
/* math.h under MacOS defines y1 and j1 as global variable */
gd->gix = gix;
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 (int i = 0; i < GeoIndexFIXEDPOINTS; i++) {
double xx1 = (gix->fixed.x)[i];
double yy1 = (gix->fixed.y)[i];
double z1 = (gix->fixed.z)[i];
double snmd = (xx1 - gd->x) * (xx1 - gd->x) + (yy1 - gd->y) * (yy1 - gd->y) +
(z1 - gd->z) * (z1 - gd->z);
(gd->fixdist)[i] = (GeoFix)(asin(sqrt(snmd) / 2.0) * ARCSINFIX);
}
}
/* =================================================== */
/* GeoMetersToSNMD */
/* When searching for a point "by distance" rather than*/
/* by count, this routine is used to reverse-engineer */
/* the distance in meters into a Squared Normalized */
/* Mole Distance (SNMD), since this is faster to */
/* compute for each individual point. Hence, rather */
/* than convert all the distances to meters and compare*/
/* the system works backwards a bit so that, for each */
/* point considered, only half of the distance */
/* calculation needs to be done. This is, of course */
/* considerably faster. */
/* =================================================== */
double GeoMetersToSNMD(double meters) {
double angle, hnmd;
angle = 0.5 * meters / EARTHRADIAN;
hnmd = sin(angle); /* half normalized mole distance */
if (angle >= M_PI / 2.0)
return 4.0;
else
return hnmd * hnmd * 4.0;
}
double GeoFixtoSNMD(GeoFix gf) {
double x;
x = gf;
x = x / ARCSINFIX;
x = sin(x);
x = x * x;
x = x * 4.0;
return x;
}
/* =================================================== */
/* GeoSetDistance */
/* During a search (of either type), the target point */
/* is first "detailed". When the distance of interest */
/* to the target point is known (either at the start */
/* of a search-by-distance or each time a new good */
/* point is found during a search-by-count) this */
/* routine is called to set the snmd and distrej valeus*/
/* so that as much as possible is known to speed up */
/* consideration of any new points */
/* =================================================== */
void GeoSetDistance(GeoDetailedPoint* gd, double snmd) {
GeoFix gf;
int i;
gd->snmd = snmd;
gf = (GeoFix)(asin(sqrt(snmd) / 2.0) * ARCSINFIX);
gf++;
for (i = 0; i < GeoIndexFIXEDPOINTS; i++) {
if ((gd->fixdist)[i] <= gf)
(gd->distrej)[i] = 0;
else
(gd->distrej)[i] = (gd->fixdist)[i] - gf;
}
}
/* =================================================== */
/* GeoStackSet routine */
/* The searches (by count and by distance) both start */
/* by detailing the point and then calling GeoStackSet */
/* Starting from the root pot (pot 1) the tree is */
/* descended towards the (actually the earliest) pot */
/* that could contain the target point. As the */
/* descent proceeds, the other child of each parent pot*/
/* is put onto the stack, so that after the routine */
/* completes, the pots on the stack are a division of */
/* the index into a set of (disjoint) intervals with */
/* a strong tendency for the ones containing near */
/* points (on the Hilbert curve, anyway) to be on the */
/* to of the stack and to contain few points */
/* =================================================== */
void GeoStackSet(GeoStack* gk, GeoDetailedPoint* gd, GeoResults* gr) {
int pot;
GeoIx* gix;
GeoPot* gp;
gix = gd->gix;
gk->gr = gr;
gk->gd = gd;
gk->stacksize = 0;
pot = 1;
while (1) {
gp = gix->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;
}
/* =================================================== */
/* GeoResultsCons routine */
/* Constructs (allocates) a new structure suitable for */
/* holding the results of a search. The GeoResults */
/* structure just holds the slotid of each point chosen*/
/* and the (SNMD) distance to the target point */
/* =================================================== */
GeoResults* GeoResultsCons(int alloc) {
GeoResults* gres;
int* sa;
double* dd;
if (alloc <= 0) {
return NULL;
}
gres = static_cast<GeoResults*>(
TRI_Allocate(TRI_UNKNOWN_MEM_ZONE, sizeof(GeoResults), false));
sa = static_cast<int*>(
TRI_Allocate(TRI_UNKNOWN_MEM_ZONE, alloc * sizeof(int), false));
dd = static_cast<double*>(
TRI_Allocate(TRI_UNKNOWN_MEM_ZONE, alloc * sizeof(double), false));
if ((gres == NULL) || (sa == NULL) || (dd == NULL)) {
if (gres != NULL) {
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gres);
}
if (sa != NULL) {
TRI_Free(TRI_UNKNOWN_MEM_ZONE, sa);
}
if (dd != NULL) {
TRI_Free(TRI_UNKNOWN_MEM_ZONE, dd);
}
return NULL;
}
gres->pointsct = 0;
gres->allocpoints = alloc;
gres->slot = sa;
gres->snmd = dd;
/* no need to initialize maxsnmd */
return gres;
}
/* =================================================== */
/* GeoResultsStartCount */
/* The GeoResultsCons routine allocates the memory */
/* but if the search is by count, it is also necessary */
/* to initialize the results list with "fake" points */
/* at the impossible SNMD distance of 10, so that any */
/* real point will be closer than that and be taken */
/* The GeoResultsStartCount routine does just that */
/* =================================================== */
void GeoResultsStartCount(GeoResults* gr) {
int i;
for (i = 0; i < gr->allocpoints; i++) {
gr->slot[i] = 0;
gr->snmd[i] = 10.0;
}
}
/* =================================================== */
/* GeoResultsInsertPoint */
/* when a point is to be considered as a candidate for */
/* being returned in a search-by-count process, the */
/* slot and snmd are presented to this routine. If the*/
/* point is too distant, it is ignored. Otherwise the */
/* most distant "old" point (which is always at zero */
/* as the results are maintained as a priority queue */
/* in this case) is discarded, and the new point must */
/* be put into its proper place to re-establish the */
/* priority queue - that every entry n is greater than */
/* or equal, in SNMD distance, than both its children */
/* which are at 2n+1 and 2n+2 */
/* =================================================== */
void GeoResultsInsertPoint(GeoResults* gr, int slot, double snmd) {
/* math.h under MacOS defines y1 and j1 as global variable */
int i, temp;
if (snmd >= gr->snmd[0]) return;
if (gr->slot[0] == 0) gr->pointsct++;
i = 0; /* i is now considered empty */
while (1) {
int jj1 = 2 * i + 1;
int jj2 = 2 * i + 2;
if (jj1 < gr->allocpoints) {
if (jj2 < gr->allocpoints) {
if (gr->snmd[jj1] > gr->snmd[jj2]) {
temp = jj1;
// jj1=jj2;
jj2 = temp;
}
/* so now jj2 is >= jj1 */
if (gr->snmd[jj2] <= snmd) {
gr->snmd[i] = snmd;
gr->slot[i] = slot;
return;
}
gr->snmd[i] = gr->snmd[jj2];
gr->slot[i] = gr->slot[jj2];
i = jj2;
continue;
}
if (gr->snmd[jj1] <= snmd) {
gr->snmd[i] = snmd;
gr->slot[i] = slot;
return;
}
gr->snmd[i] = gr->snmd[jj1];
gr->slot[i] = gr->slot[jj1];
i = jj1;
continue;
}
gr->snmd[i] = snmd;
gr->slot[i] = slot;
return;
}
}
/* =================================================== */
/* GeoResultsGrow */
/* During a search-by distance (the search-by-count */
/* allocates the correct size at the outset) it may be */
/* necessary to return an unbounded amount of data. */
/* initially 100 entries are allocted, but this routine*/
/* ensures that another one is available. If the */
/* allocation fails, -1 is returned. */
/* =================================================== */
int GeoResultsGrow(GeoResults* gr) {
int newsiz;
int* sa;
double* dd;
if (gr->pointsct < gr->allocpoints) return 0;
/* otherwise grow by about 50% */
newsiz = gr->pointsct + (gr->pointsct / 2) + 1;
if (newsiz > 1000000000) return -1;
sa = static_cast<int*>(
TRI_Reallocate(TRI_UNKNOWN_MEM_ZONE, gr->slot, newsiz * sizeof(int)));
dd = static_cast<double*>(
TRI_Reallocate(TRI_UNKNOWN_MEM_ZONE, gr->snmd, newsiz * sizeof(double)));
if ((sa == NULL) || (dd == NULL)) {
if (sa != NULL) gr->slot = sa;
if (dd != NULL) gr->snmd = dd;
return -1;
}
gr->slot = sa;
gr->snmd = dd;
gr->allocpoints = newsiz;
return 0;
}
/* =================================================== */
/* GeoAnswers */
/* At the end of any search (of either type) the */
/* GeoResults structure holds the slotid and snmd */
/* distance of the points to be returned. This routine*/
/* constructs and populates the GeoCoordinates */
/* structure with the require data by fetching the */
/* coodinates from the index, and by convertin the */
/* snmd distance into meters. It should be noticed */
/* that the latitude and longitude are copied into the */
/* new data, so that the GeoCoordinates structure */
/* remains valid even if the index is subsequently */
/* updated or even freed. NOTICE also that the */
/* distances returned may not agree precisely with the */
/* distances that could be calculated by a separate */
/* call to GeoIndex_distance because of rounding errors*/
/* =================================================== */
GeoCoordinates* GeoAnswers(GeoIx* gix, GeoResults* gr) {
GeoCoordinates* ans;
GeoCoordinate* gc;
int i, j;
double mole;
if (gr->pointsct == 0) {
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gr->slot);
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gr->snmd);
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gr);
return NULL;
}
ans = static_cast<GeoCoordinates*>(
TRI_Allocate(TRI_UNKNOWN_MEM_ZONE, sizeof(GeoCoordinates), false));
gc = static_cast<GeoCoordinate*>(TRI_Allocate(
TRI_UNKNOWN_MEM_ZONE, gr->pointsct * sizeof(GeoCoordinate), false));
if ((ans == NULL) || (gc == NULL)) {
if (ans != NULL) {
TRI_Free(TRI_UNKNOWN_MEM_ZONE, ans);
}
if (gc != NULL) {
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gc);
}
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gr->slot);
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gr->snmd);
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gr);
return NULL;
}
ans->length = gr->pointsct;
ans->coordinates = gc;
j = 0;
for (i = 0; i < gr->allocpoints; i++) {
if (j >= gr->pointsct) break;
int slot = gr->slot[i];
if (slot == 0) continue;
ans->coordinates[j].latitude = (gix->gc)[slot].latitude;
ans->coordinates[j].longitude = (gix->gc)[slot].longitude;
ans->coordinates[j].data = (gix->gc)[slot].data;
mole = sqrt(gr->snmd[i]);
if (mole > 2.0) mole = 2.0; /* make sure arcsin succeeds! */
gr->snmd[j] = 2.0 * EARTHRADIAN * asin(mole / 2.0);
j++;
}
ans->distances = gr->snmd;
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gr->slot);
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gr);
return ans;
}
/* =================================================== */
/* GeoPotJunk */
/* A detailed point containing the target point set */
/* with the current distance is compared to a pot */
/* If any of the fixed points are too close to all the */
/* descendents of a pot, 1 is returned to indicate that*/
/* the pot is "junk" = it may be ignored in its */
/* entirety because it contains no points close enough */
/* to the target. Otherwise 0 is returned. */
/* =================================================== */
int GeoPotJunk(GeoDetailedPoint* gd, int pot) {
int i;
GeoPot* gp;
gp = (gd->gix)->pots + pot;
for (i = 0; i < GeoIndexFIXEDPOINTS; i++)
if (gp->maxdist[i] < gd->distrej[i]) return 1;
return 0;
}
/* =================================================== */
/* GeoSNMD */
/* Finds the SNMD (Squared NormalizedMole Distance) */
/* from the point (which must be "detailed" gd, and the*/
/* ordinary point (just given by lat/longitude) */
/* The cartesian coordinates of the ordinary point are */
/* found, and then the differences squared returned. */
/* =================================================== */
double GeoSNMD(GeoDetailedPoint* gd, GeoCoordinate* c) {
double x, y, z;
z = sin(c->latitude * 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);
}
/* =================================================== */
/* GeoIndex_PointsWithinRadius */
/* This is the basic user-visible call to find all the */
/* the points in the index that are within the */
/* specified distance of the target point */
/* First the GeoIndex must be cast to the correct */
/* (GeoIx) structure so that it can be used! */
/* the result structure is then set up initially to */
/* hold up to 100 results points, and the point is then*/
/* detailed (GeoString, x,y,z and distances to fixed */
/* points). The stack is then populated with the */
/* initial descending set of pots ending with the one */
/* nearest the target point, and the distance set on */
/* the detailed point by converting the meters into an */
/* SNMD. The pots on the stack are then considered. */
/* If the call to GeoPotJunk indicates that there are */
/* no points in that pot within the required circle, */
/* the pot is discarded. Otherwise, if the pot is a */
/* leaf pot, the points are considered individually, */
/* and notice the recovery to free everything if there */
/* is a need to grow the results structure and there */
/* is not enough memory. If the pot is not a leaf pot */
/* it is replaced on the stack by both its children */
/* Processing continues until the stack is empty */
/* At the end, the GeoAnswers routine is used to */
/* convert the pot/snmd collection of the GeoResults */
/* structure, into the distance (in meters) and the */
/* GeoCoordinate data (lat/longitude and data pointer) */
/* needed for the return to the caller. */
/* =================================================== */
GeoCoordinates* GeoIndex_PointsWithinRadius(GeoIdx* gi, GeoCoordinate* c,
double d) {
GeoResults* gres;
GeoCoordinates* answer;
GeoDetailedPoint gd;
GeoStack gk;
GeoPot* gp;
int r, slot, i;
double snmd, maxsnmd;
GeoIx* gix;
if (c->longitude < -180.0) return NULL;
if (c->longitude > 180.0) return NULL;
if (c->latitude < -90.0) return NULL;
if (c->latitude > 90.0) return NULL;
gix = (GeoIx*)gi;
gres = GeoResultsCons(100);
if (gres == NULL) return NULL;
GeoMkDetail(gix, &gd, c);
GeoStackSet(&gk, &gd, gres);
maxsnmd = GeoMetersToSNMD(d);
GeoSetDistance(&gd, maxsnmd);
gk.stacksize++;
while (gk.stacksize >= 1) {
gk.stacksize--;
int pot = gk.potid[gk.stacksize];
if (GeoPotJunk(&gd, pot)) continue;
gp = gix->pots + pot;
if (gp->LorLeaf == 0) {
for (i = 0; i < gp->RorPoints; i++) {
slot = gp->points[i];
snmd = GeoSNMD(&gd, gix->gc + slot);
if (snmd > (maxsnmd * 1.00000000000001)) continue;
r = GeoResultsGrow(gres);
if (r == -1) {
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gres->snmd);
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gres->slot);
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gres);
return NULL;
}
gres->slot[gres->pointsct] = slot;
gres->snmd[gres->pointsct] = snmd;
gres->pointsct++;
}
} else {
gk.potid[gk.stacksize++] = gp->LorLeaf;
gk.potid[gk.stacksize++] = gp->RorPoints;
}
}
answer = GeoAnswers(gix, gres);
return answer; /* note - this may be NULL */
}
/* =================================================== */
/* GeoIndex_NearestCountPoints */
/* The other user-visible search call, which finds the */
/* nearest <count> points for a user-specified <count> */
/* processing is not dissimilar to the previous routine*/
/* but here the results structure is allocated at the */
/* correct size and used as a priority queue. Since */
/* it always helps if more points are found (the */
/* distance of interest drops, so that pots are more */
/* readily rejected) some care is taken when a pot is */
/* not rejected to put the one most likely to contain */
/* useful points onto the top of the stack for early */
/* processing. */
/* =================================================== */
GeoCoordinates* GeoIndex_NearestCountPoints(GeoIdx* gi, GeoCoordinate* c,
int count) {
GeoResults* gr;
GeoDetailedPoint gd;
GeoCoordinates* answer;
GeoStack gk;
GeoPot* gp;
int slot, i, left;
double snmd;
GeoIx* gix;
if (c->longitude < -180.0) return NULL;
if (c->longitude > 180.0) return NULL;
if (c->latitude < -90.0) return NULL;
if (c->latitude > 90.0) return NULL;
gix = (GeoIx*)gi;
gr = GeoResultsCons(count);
if (gr == NULL) return NULL;
GeoMkDetail(gix, &gd, c);
GeoStackSet(&gk, &gd, gr);
GeoResultsStartCount(gr);
left = count;
while (gk.stacksize >= 0) {
int pot = gk.potid[gk.stacksize--];
gp = gix->pots + pot;
if (left <= 0) {
GeoSetDistance(&gd, gr->snmd[0]);
if (GeoPotJunk(&gd, pot)) continue;
}
if (gp->LorLeaf == 0) {
for (i = 0; i < gp->RorPoints; i++) {
slot = gp->points[i];
snmd = GeoSNMD(&gd, gix->gc + slot);
GeoResultsInsertPoint(gr, slot, snmd);
left--;
if (left < -1) left = -1;
}
} else {
if (gd.gs > gp->middle) {
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(gix, gr);
return answer; /* note - this may be NULL */
}
/* =================================================== */
/* GeoIndexFreeSlot */
/* return the specified slot to the free list */
/* =================================================== */
void GeoIndexFreeSlot(GeoIx* gix, int slot) {
gix->gc[slot].latitude = gix->gc[0].latitude;
gix->gc[0].latitude = slot;
}
/* =================================================== */
/* GeoIndexNewSlot */
/* If there is a fre slot already on the free list, */
/* just return its slot number. Otherwise the entire */
/* slot list is realloc'd. Although this might change */
/* the physical memory location of all the indexed */
/* points, this is not a problem since the slotid */
/* values are not changed. */
/* The GeoIndexGROW, which specifies the percentage */
/* of growth to be used, is in GeoIndex.h. Notice also*/
/* that some care is take to ensure that, in the case */
/* of memory allocation failure, the index is still */
/* kept unchanged even though the new point cannot be */
/* added to the index. */
/* =================================================== */
int GeoIndexNewSlot(GeoIx* gix) {
int j;
GeoCoordinate* gc;
if (gix->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 */
long long x = gix->slotct;
long long y = 100 + GeoIndexGROW;
x = x * y + 99;
y = 100;
x = x / y;
if (x > 2000000000L) return -2;
int newslotct = (int)x;
gc = static_cast<GeoCoordinate*>(TRI_Reallocate(
TRI_UNKNOWN_MEM_ZONE, gix->gc, newslotct * sizeof(GeoCoordinate)));
if (gc == NULL) {
return -2;
}
gix->gc = gc;
// update memory usage
gix->_memoryUsed -= gix->slotct * sizeof(GeoCoordinate);
gix->_memoryUsed += newslotct * sizeof(GeoCoordinate);
for (j = gix->slotct; j < newslotct; j++) GeoIndexFreeSlot(gix, j);
gix->slotct = newslotct;
}
j = (int)(gix->gc[0].latitude);
gix->gc[0].latitude = gix->gc[j].latitude;
return j;
}
/* =================================================== */
/* GeoFind */
/* This routine is used during insertion and removal, */
/* but is not used during the searches. */
/* Find the given point if it is in the index, and set */
/* the GeoPath data structure to give the path from the*/
/* root pot (pot 1) to the leaf pot, if any, containing*/
/* the sepecified (detailed) point, or - if the point */
/* is not present, to the first leaf pot into which the*/
/* specified point may be inserted. */
/* To start with, the index tree is descended, starting*/
/* with the root (which, rather bizzarly, is at the */
/* top of this tree!) always taking the right branch if*/
/* both would do, to reach the rightmost leaf pot that */
/* could contain the specified point. */
/* We then proceed leftwards through the points until */
/* either the specified point is found in the index, or*/
/* the first leaf pot is found that could contain the */
/* specified point. It is worth noting that the first */
/* pot of all has "low-values" as its "start" GeoString*/
/* so that this process cannot go off the front of the */
/* index. Notice also that it is not expected to be */
/* very common that a large number of points with the */
/* same GeoString (so within 30 centimeters!) will be */
/* inserted into the index, and that even if there are */
/* the inefficiency of this code is only moderate, and */
/* manifests itself only during maintenance */
/* the return value is 1 if the point is found and 2 */
/* if it is not found */
/* =================================================== */
int GeoFind(GeoPath* gt, GeoDetailedPoint* gd) {
int pot, pot1;
int i;
int slot;
GeoIx* gix;
GeoCoordinate* gc;
GeoPot* gp;
gix = gd->gix;
gt->gix = gix;
pot = 1;
gt->pathlength = 0;
while (1) {
gp = gix->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; i < gp->RorPoints; i++) {
slot = gp->points[i];
gc = gix->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->start < gd->gs) 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 = gix->pots + pot1;
if (pot == gp->RorPoints) break; /* cannot go off the front */
pot = pot1;
}
gp = gix->pots + pot1;
pot = gp->LorLeaf;
/* now we have a pot whose iterated right child we want */
while (1) {
gp = gix->pots + pot;
gt->path[gt->pathlength] = pot;
gt->pathlength++;
if (gp->LorLeaf == 0) break;
pot = gp->RorPoints;
}
}
return 2;
}
/* =================================================== */
/* GeoPopulateMaxdist */
/* During maintencance, when the points in a leaf pot */
/* have been changed, this routine merely looks at all */
/* the points in the pot, details them, and rebuilds */
/* the list of maximum distances. */
/* =================================================== */
void GeoPopulateMaxdist(GeoIx* gix, GeoPot* gp, GeoString* gsa) {
int i, j;
GeoDetailedPoint gd;
gsa[0] = 0x1FFFFFFFFFFFFFll;
gsa[1] = 0ll;
for (j = 0; j < GeoIndexFIXEDPOINTS; j++) gp->maxdist[j] = 0;
for (i = 0; i < gp->RorPoints; i++) {
GeoMkDetail(gix, &gd, gix->gc + gp->points[i]);
for (j = 0; j < GeoIndexFIXEDPOINTS; j++)
if (gd.fixdist[j] > gp->maxdist[j]) gp->maxdist[j] = gd.fixdist[j];
if (gd.gs < gsa[0]) gsa[0] = gd.gs;
if (gd.gs > gsa[1]) gsa[1] = gd.gs;
}
gp->level = 1;
}
/* =================================================== */
/* GeoGetPot */
/* This routine simply converts a path and a height */
/* into a pot id. */
/* =================================================== */
int GeoGetPot(GeoPath* gt, int height) {
return gt->path[gt->pathlength - height];
}
/* =================================================== */
/* GeoAdjust */
/* During insertion and deletion, this routine is used */
/* to populate the data correctly for the parent pot */
/* specified (which may not be a leaf pot) by taking */
/* the data from the child pots. It populates the */
/* start, middle and end GeoStrings, the level, and */
/* the maximum distances to the fixed points. */
/* =================================================== */
void GeoAdjust(GeoIx* gix, int potx) /* the kids are alright */
{
int poty, potz; /* x = (yz) */
int i;
GeoPot* gpx;
GeoPot* gpy;
GeoPot* gpz;
gpx = gix->pots + potx;
poty = gpx->LorLeaf;
gpy = gix->pots + poty;
potz = gpx->RorPoints;
gpz = gix->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; i < GeoIndexFIXEDPOINTS; i++) {
gpx->maxdist[i] = gpy->maxdist[i];
if (gpx->maxdist[i] < gpz->maxdist[i]) gpx->maxdist[i] = gpz->maxdist[i];
}
}
/* =================================================== */
/* RotateLeft */
/* The operation used during tree balancing to convert */
/* A(BC) into (AB)C. To start with, E is A(BC) and */
/* D is BC. D is then change to be (AB) and */
/* GeoAdjust is used to re-populate its data. E is */
/* then set to be DC = (AB)C, and again GeoAdjust is */
/* used to set the GeoStrings, level and distances to */
/* the fixed points, taking the data from the children */
/* in both cases */
/* =================================================== */
void RotateLeft(GeoIx* gix, int pote) {
int pota, potb, potc, potd;
GeoPot* gpd;
GeoPot* gpe;
gpe = gix->pots + pote;
potd = gpe->RorPoints;
gpd = gix->pots + potd;
pota = gpe->LorLeaf;
potb = gpd->LorLeaf;
potc = gpd->RorPoints;
gpd->LorLeaf = pota;
gpd->RorPoints = potb;
GeoAdjust(gix, potd);
gpe->LorLeaf = potd;
gpe->RorPoints = potc;
GeoAdjust(gix, pote);
}
/* =================================================== */
/* RotateRight */
/* The mirror-image or inverse of RotateLeft. */
/* Changes (AB)C into A(BC). The given parent pot is */
/* E = (AB)C and D is AB. D is then reused to be BC */
/* and GeoAdjusted, and then E set to be AD = A(BC) and*/
/* also GeoAdjusted */
/* =================================================== */
void RotateRight(GeoIx* gix, int pote) {
int pota, potb, potc, potd;
GeoPot* gpd;
GeoPot* gpe;
gpe = gix->pots + pote;
potd = gpe->LorLeaf;
gpd = gix->pots + potd;
pota = gpd->LorLeaf;
potb = gpd->RorPoints;
potc = gpe->RorPoints;
gpd->LorLeaf = potb;
gpd->RorPoints = potc;
GeoAdjust(gix, potd);
gpe->LorLeaf = pota;
gpe->RorPoints = potd;
GeoAdjust(gix, pote);
}
/* =================================================== */
/* GeoIndex_insert */
/* The user-facing routine to insert a new point into */
/* the index. First the index is cast into a GeoIx */
/* so that it can be used, and then the point is */
/* sanity checked. The point is then detailed and the */
/* GeoFind routine called. If the point is found, this*/
/* is an error. Otherwise a new slot is populated with*/
/* the data from the point, and then the point is put */
/* into the first leaf pot into which it may go based */
/* on its GeoString value. If there is no room in that*/
/* pot, the pot is split into two (necessitating a tree*/
/* balancing operation) which starts by obtaining the */
/* two new pots. . . continued below */
/* =================================================== */
int GeoIndex_insert(GeoIdx* gi, GeoCoordinate* c) {
int i, j, slot, pot, pot1;
int pota, poty, potz;
int 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];
GeoIx* gix;
gix = (GeoIx*)gi;
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(gix, &gd, c);
i = GeoFind(&gt, &gd);
if (i == 1) return -1;
pot = gt.path[gt.pathlength - 1];
gp = gix->pots + pot;
/* new point, so we try to put it in */
slot = GeoIndexNewSlot(gix);
if (slot == -2) return -2; /* no room :( */
gix->gc[slot].latitude = c->latitude;
gix->gc[slot].longitude = c->longitude;
gix->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(gix);
int pot2 = GeoIndexNewPot(gix);
gp = gix->pots + pot; /* may have re-alloced! */
if ((pot1 == -2) || (pot2 == -2)) {
GeoIndexFreeSlot(gix, slot);
if (pot1 != -2) GeoIndexFreePot(gix, pot1);
if (pot2 != -2) GeoIndexFreePot(gix, pot2);
return -2;
}
/* =================================================== */
/* GeoIndex_insert continued */
/* New pots are pot1 and pot2 which will be the new */
/* leaf pots with half the points each, and the old */
/* pot will become the parent of both of them */
/* After moving all the points to pot2, the half with */
/* the lowest GeoString are moved into pot1. The two */
/* pots are then inspected with GeoPopulateMaxdist */
/* to ascertain what the actual distances and GeoString*/
/* values are. The GeoString boundary between the two */
/* pots is set at the midpoint between the current */
/* actual boundaries and finally the current pot is */
/* set to be either pot1 or pot2 depending on where the*/
/* new point (which has still not been inserted) shoud */
/* go. Continued below . . . . */
/* =================================================== */
gp1 = gix->pots + pot1;
gp2 = gix->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; i < gp->RorPoints; i++) gp2->points[i] = gp->points[i];
/* move the first half of the points from pot2 to pot1 */
GeoString mings;
for (i = 0; i < (GeoIndexPOTSIZE / 2); i++) {
mings = 0x1FFFFFFFFFFFFFll;
int js = 0;
for (j = 0; j < gp2->RorPoints; j++) {
GeoString gs = GeoMkHilbert(gix->gc + gp2->points[j]);
if (gs < mings) {
mings = gs;
js = j;
}
}
gp1->points[gp1->RorPoints] = gp2->points[js];
gp2->points[js] = gp2->points[gp2->RorPoints - 1];
gp2->RorPoints--;
gp1->RorPoints++;
}
GeoPopulateMaxdist(gix, gp2, gsa);
mings = gsa[0];
GeoPopulateMaxdist(gix, 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(gix, pot);
gt.pathlength++;
if (gd.gs < mings) {
gp = gp1;
gt.path[gt.pathlength - 1] = pot1;
} else {
gp = gp2;
gt.path[gt.pathlength - 1] = pot2;
}
}
/* =================================================== */
/* GeoIndex_insert continued */
/* finally the new point is inserted into the pot, and */
/* the maximum distances to the fixed points propogated*/
/* up as far as necessary. The rebalancing of the tree*/
/* is then done, but only if the pot splitting happend */
/* to rebalance, the sequence of pots going back up is */
/* traversed using the path structure, and the standard*/
/* AVL balancing is used by doing the necessary */
/* rotations and level changes necessary to ensure that*/
/* every parent has at least one child one level lower */
/* and the other child is either also one level lower, */
/* or two levels lower. The details are also given in */
/* the accompanying documentation */
/* =================================================== */
/* so we have a pot and a path we can use */
/* gp is the pot, gt set correctly */
gp->points[gp->RorPoints] = slot;
gp->RorPoints++;
/* now propagate the maxdistances */
for (i = 0; i < GeoIndexFIXEDPOINTS; i++) {
j = gt.pathlength - 1;
while (j >= 0) {
if (gd.fixdist[i] > gix->pots[gt.path[j]].maxdist[i])
gix->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 (true) {
int potx = GeoGetPot(&gt, height);
gpx = gix->pots + potx;
int lvx = gpx->level;
if (potx == 1) break;
/* root pot ? */
pot1 = GeoGetPot(&gt, height + 1); /* pot1=parent(x) */
gp1 = gix->pots + pot1;
int lv1 = gp1->level;
if (lv1 > lvx) break;
if (gp1->LorLeaf == potx) /* gpx is the left child? */
{
pota = gp1->RorPoints; /* 1 = (xa) */
gpa = gix->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 = gix->pots + poty;
lvy = gpy->level;
potz = gpx->LorLeaf;
gpz = gix->pots + potz;
lvz = gpz->level;
if (lvy <= lvz) {
RotateRight(gix, pot1);
height++;
continue;
}
RotateLeft(gix, potx);
RotateRight(gix, pot1);
} else /* gpx is the right child */
{
pota = gp1->LorLeaf; /* 1 = (ax) */
gpa = gix->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 = gix->pots + poty;
lvy = gpy->level;
potz = gpx->RorPoints;
gpz = gix->pots + potz;
lvz = gpz->level;
if (lvy <= lvz) {
RotateLeft(gix, pot1);
height++;
continue;
}
RotateRight(gix, potx);
RotateLeft(gix, pot1);
}
}
return 0;
}
/* =================================================== */
/* GeoIndex_remove */
/* As a user-facing routine, this starts by casting the*/
/* GeoIndex structure to a GeoIx, so that its members */
/* can be accessed. The point is then detailed, and */
/* GeoFind is used to check whether it is there. If */
/* not, this is an error. Otherwise the point is */
/* removed from the pot and the distances recalculated */
/* using the GeoPopulateMaxdist routine. It is then */
/* checked whether there are now too few points in the */
/* pot that used to contain the point, and if so there */
/* are eight cases as to what is to be done. In four */
/* of them, a point is moved from the adjacent leaf pot*/
/* which may be at the same level or one lower, and may*/
/* be either side of the current one. This is done if */
/* there are too many points in the two leaf pots to */
/* amalgamate them. In the other four cases the two */
/* leaf pots are amalgamated, which results in the */
/* releasing of two pots (which are put back into the */
/* free chain using GeoIndexFreePot) Continued . . . . */
/* =================================================== */
int GeoIndex_remove(GeoIdx* gi, GeoCoordinate* c) {
GeoDetailedPoint gd;
int rebalance;
int levn, levc;
GeoPot* gp;
int potp;
GeoPot* gpp;
int potb;
GeoPot* gpb;
int potn;
GeoPot* gpn;
int potc;
GeoPot* gpc;
GeoPath gt;
GeoString gsa[2];
int i, pot, potix, slot, pathix;
GeoIx* gix;
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;
gix = (GeoIx*)gi;
GeoMkDetail(gix, &gd, c);
i = GeoFind(&gt, &gd);
if (i != 1) return -1;
pot = gt.path[gt.pathlength - 1];
gp = gix->pots + pot;
potix = gt.path[gt.pathlength];
slot = gp->points[potix];
GeoIndexFreeSlot(gix, slot);
gp->points[potix] = gp->points[gp->RorPoints - 1];
gp->RorPoints--;
GeoPopulateMaxdist(gix, gp, gsa);
if (pot == 1) return 0; /* just allow root pot to have fewer points */
rebalance = 0;
if ((2 * gp->RorPoints) < GeoIndexPOTSIZE) {
int j, js;
GeoString mings, gs;
potp = gt.path[gt.pathlength - 2];
gpp = gix->pots + potp;
if (gpp->LorLeaf == pot) {
/* Left */
potb = gpp->RorPoints;
gpb = gix->pots + potb;
if (gpb->LorLeaf == 0) {
/* Left Brother */
if ((gpb->RorPoints + gp->RorPoints) > GeoIndexPOTSIZE) {
/* Left Brother Lots */
mings = 0x1FFFFFFFFFFFFFll;
js = 0;
for (j = 0; j < gpb->RorPoints; j++) {
gs = GeoMkHilbert(gix->gc + gpb->points[j]);
if (gs < mings) {
mings = gs;
js = j;
}
}
gp->points[gp->RorPoints] = gpb->points[js];
gpb->points[js] = gpb->points[gpb->RorPoints - 1];
gpb->RorPoints--;
gp->RorPoints++;
GeoPopulateMaxdist(gix, gp, gsa);
mings = gsa[1];
GeoPopulateMaxdist(gix, gpb, gsa);
mings = (mings + gsa[0]) / 2ll;
gp->end = mings;
gpb->start = mings;
gpp->middle = mings;
GeoAdjust(gix, potp);
} else {
/* Left Brother Few */
gpp->LorLeaf = 0;
i = 0;
for (j = 0; j < gpb->RorPoints; j++)
gpp->points[i++] = gpb->points[j];
for (j = 0; j < gp->RorPoints; j++) gpp->points[i++] = gp->points[j];
gpp->RorPoints = i;
GeoIndexFreePot(gix, pot);
GeoIndexFreePot(gix, potb);
GeoPopulateMaxdist(gix, gpp, gsa);
gt.pathlength--;
rebalance = 1;
}
} else {
/* Left Nephew */
potn = gpb->LorLeaf;
gpn = gix->pots + potn;
if ((gpn->RorPoints + gp->RorPoints) > GeoIndexPOTSIZE) {
/* Left Nephew Lots */
mings = 0x1FFFFFFFFFFFFFll;
js = 0;
for (j = 0; j < gpn->RorPoints; j++) {
gs = GeoMkHilbert(gix->gc + gpn->points[j]);
if (gs < mings) {
mings = gs;
js = j;
}
}
gp->points[gp->RorPoints] = gpn->points[js];
gpn->points[js] = gpn->points[gpn->RorPoints - 1];
gpn->RorPoints--;
gp->RorPoints++;
GeoPopulateMaxdist(gix, gp, gsa);
mings = gsa[1];
GeoPopulateMaxdist(gix, gpn, gsa);
mings = (mings + gsa[0]) / 2ll;
gp->end = mings;
gpn->start = mings;
gpb->start = mings;
gpp->middle = mings;
GeoAdjust(gix, potb);
GeoAdjust(gix, potp);
} else {
/* Left Nephew Few */
potc = gpb->RorPoints;
i = gp->RorPoints;
for (j = 0; j < gpn->RorPoints; j++) gp->points[i++] = gpn->points[j];
gp->RorPoints = i;
gpp->RorPoints = potc;
gpp->middle = gpb->middle;
gp->end = gpp->middle;
GeoIndexFreePot(gix, potn);
GeoIndexFreePot(gix, potb);
GeoPopulateMaxdist(gix, gp, gsa);
GeoAdjust(gix, potp);
gt.pathlength--;
rebalance = 1;
}
}
} else {
/* Right */
potb = gpp->LorLeaf;
gpb = gix->pots + potb;
if (gpb->LorLeaf == 0) {
/* Right Brother */
if ((gpb->RorPoints + gp->RorPoints) > GeoIndexPOTSIZE) {
/* Right Brother Lots */
mings = 0ll;
js = 0;
for (j = 0; j < gpb->RorPoints; j++) {
gs = GeoMkHilbert(gix->gc + gpb->points[j]);
if (gs > mings) {
mings = gs;
js = j;
}
}
gp->points[gp->RorPoints] = gpb->points[js];
gpb->points[js] = gpb->points[gpb->RorPoints - 1];
gpb->RorPoints--;
gp->RorPoints++;
GeoPopulateMaxdist(gix, gp, gsa);
mings = gsa[0];
GeoPopulateMaxdist(gix, gpb, gsa);
mings = (mings + gsa[1]) / 2ll;
gp->start = mings;
gpb->end = mings;
gpp->middle = mings;
GeoAdjust(gix, potp);
} else {
/* Right Brother Few */
/* observe this is identical to Left Brother Few */
gpp->LorLeaf = 0;
i = 0;
for (j = 0; j < gpb->RorPoints; j++)
gpp->points[i++] = gpb->points[j];
for (j = 0; j < gp->RorPoints; j++) gpp->points[i++] = gp->points[j];
gpp->RorPoints = i;
GeoIndexFreePot(gix, pot);
GeoIndexFreePot(gix, potb);
GeoPopulateMaxdist(gix, gpp, gsa);
gt.pathlength--;
rebalance = 1;
}
} else {
/* Right Nephew */
potn = gpb->RorPoints;
gpn = gix->pots + potn;
if ((gpn->RorPoints + gp->RorPoints) > GeoIndexPOTSIZE) {
/* Right Nephew Lots */
mings = 0ll;
js = 0;
for (j = 0; j < gpn->RorPoints; j++) {
gs = GeoMkHilbert(gix->gc + gpn->points[j]);
if (gs > mings) {
mings = gs;
js = j;
}
}
gp->points[gp->RorPoints] = gpn->points[js];
gpn->points[js] = gpn->points[gpn->RorPoints - 1];
gpn->RorPoints--;
gp->RorPoints++;
GeoPopulateMaxdist(gix, gp, gsa);
mings = gsa[0];
GeoPopulateMaxdist(gix, gpn, gsa);
mings = (mings + gsa[1]) / 2ll;
gp->start = mings;
gpn->end = mings;
gpb->end = mings;
gpp->middle = mings;
GeoAdjust(gix, potb);
GeoAdjust(gix, potp);
} else {
/* Right Nephew Few */
potc = gpb->LorLeaf;
i = gp->RorPoints;
for (j = 0; j < gpn->RorPoints; j++) gp->points[i++] = gpn->points[j];
gp->RorPoints = i;
gpp->LorLeaf = potc;
gpp->middle = gpb->middle;
gp->start = gpb->middle;
GeoIndexFreePot(gix, potn);
GeoIndexFreePot(gix, potb);
GeoPopulateMaxdist(gix, gp, gsa);
GeoAdjust(gix, potp);
gt.pathlength--;
rebalance = 1;
}
}
}
}
/* =================================================== */
/* GeoIndex_remove continued */
/* Again the balancing of the tree is fairly standard */
/* and documented in the associated documentation to */
/* this routine. At every stage in this process the */
/* parent potp of the current pot may not be balanced */
/* as pot has just had its level reduced. To tell what*/
/* to do, the product i of the level differences is */
/* calculated. This should be 1 or 2, but may be 3 or */
/* 4, and in each case some further investigation soon */
/* shows what rotations and further upward balancing */
/* may be needed. continued . . . */
/* =================================================== */
pathix = gt.pathlength - 1;
while ((pathix > 0) && (rebalance == 1)) {
/* Deletion rebalancing */
rebalance = 0;
pathix--;
potp = gt.path[pathix];
gpp = gix->pots + potp;
int levp = gpp->level;
pot = gpp->LorLeaf;
potb = gpp->RorPoints;
gp = gix->pots + pot;
gpb = gix->pots + potb;
int lev = gp->level;
int levb = gpb->level;
i = (levp - lev) * (levp - levb);
if (i == 4) {
gpp->level--;
rebalance = 1;
}
if (i == 3) {
if ((levp - lev) == 3) {
potn = gpb->LorLeaf;
gpn = gix->pots + potn;
potc = gpb->RorPoints;
gpc = gix->pots + potc;
levn = gpn->level;
levc = gpc->level;
if (levn <= levc) {
RotateLeft(gix, potp);
if (levn < levc) rebalance = 1;
} else {
RotateRight(gix, potb);
RotateLeft(gix, potp);
rebalance = 1;
}
} else {
potn = gp->LorLeaf;
gpn = gix->pots + potn;
potc = gp->RorPoints;
gpc = gix->pots + potc;
levn = gpn->level;
levc = gpc->level;
if (levn >= levc) {
RotateRight(gix, potp);
if (levn > levc) rebalance = 1;
} else {
RotateLeft(gix, pot);
RotateRight(gix, potp);
rebalance = 1;
}
}
}
GeoAdjust(gix, potp);
}
/* =================================================== */
/* GeoIndex_remove continued */
/* In the case of deletion, it is not so easy to see */
/* what the new maximum distances are given the point */
/* deleted, so the GeoAdjust routine is used all the */
/* way up. */
/* =================================================== */
while (pathix > 0) {
pathix--;
pot = gt.path[pathix];
GeoAdjust(gix, pot);
}
return 0;
}
/* =================================================== */
/* GeoIndex_CoordinatesFree */
/* The user-facing routine that must be called by the */
/* user when the results of a search are finished with */
/* =================================================== */
void GeoIndex_CoordinatesFree(GeoCoordinates* clist) {
if (clist == nullptr) {
return;
}
TRI_Free(TRI_UNKNOWN_MEM_ZONE, clist->coordinates);
TRI_Free(TRI_UNKNOWN_MEM_ZONE, clist->distances);
TRI_Free(TRI_UNKNOWN_MEM_ZONE, clist);
}
/* =================================================== */
/* GeoIndex_hint does nothing! */
/* it is here for possible future compatibilty */
/* =================================================== */
int GeoIndex_hint(GeoIdx* gi, int hint) { return 0; }
/* =================================================== */
/* GeoCr structure */
/* This is the REAL GeoCursor structure - the one in */
/* the GeoIndex.h file is just a sham (it says it is */
/* a char!) to keep the structure private so that the */
/* GeoIndex.h is short and contains only data of */
/* interest to the user. */
/* =================================================== */
typedef struct {
int pot;
GeoFix dist;
} hpot; // pot for putting on the heap
bool hpotcompare(hpot a, hpot b) { return (a.dist > b.dist); }
typedef struct {
int slot;
double snmd;
} hslot; // pot for putting on the heap
bool hslotcompare(hslot a, hslot b) {
if (a.snmd > b.snmd)
return true;
else
return false;
}
typedef struct {
GeoIx* Ix; /* GeoIndex */
GeoDetailedPoint gd;
double potsnmd;
double slotsnmd;
std::vector<hpot> potheap;
std::vector<hslot> slotheap;
} GeoCr;
GeoFix makedist(GeoPot* pot, GeoDetailedPoint* gd) {
GeoFix dist, d1;
int i;
dist = 0;
for (i = 0; i < GeoIndexFIXEDPOINTS; i++) {
if (gd->fixdist[i] > pot->maxdist[i])
d1 = gd->fixdist[i] - pot->maxdist[i];
else
d1 = 0;
if (d1 > dist) dist = d1;
}
return dist;
}
GeoCursor* GeoIndex_NewCursor(GeoIdx* gi, GeoCoordinate* c) {
GeoIx* gix;
hpot hp;
if (c->longitude < -180.0) return nullptr;
if (c->longitude > 180.0) return nullptr;
if (c->latitude < -90.0) return nullptr;
if (c->latitude > 90.0) return nullptr;
gix = (GeoIx*)gi;
GeoCr* gcr = nullptr;
try {
gcr = new GeoCr;
} catch (...) {
}
if (gcr == nullptr) {
return (GeoCursor*)gcr;
}
gcr->Ix = gix;
std::vector<hpot>* p = new (&gcr->potheap) std::vector<hpot>();
std::vector<hslot>* q = new (&gcr->slotheap) std::vector<hslot>();
(void)p; // avoid compiler warnings - I just want to call
(void)q; // the constructors and have no use for p,q.
GeoMkDetail(gix, &(gcr->gd), c);
hp.pot = 1;
hp.dist = makedist(gix->pots + 1, &(gcr->gd));
gcr->potsnmd = GeoFixtoSNMD(hp.dist);
gcr->slotsnmd = 20.0;
gcr->potheap.push_back(hp);
std::push_heap(gcr->potheap.begin(), gcr->potheap.end(), hpotcompare);
return (GeoCursor*)gcr;
}
GeoCoordinates* GeoIndex_ReadCursor(GeoCursor* gc, int count) {
int i, j, r;
GeoCoordinate* ct;
GeoResults* gr;
GeoCoordinates* gcts;
GeoCr* gcr;
GeoPot pot;
int slox;
double tsnmd;
hslot hs;
hpot hp;
gcr = (GeoCr*)gc;
gr = GeoResultsCons(count);
if (gr == NULL) return NULL;
while (gr->pointsct < count) {
if (gcr->potsnmd < gcr->slotsnmd * 1.000001) {
// smash top pot - if there is one
if (gcr->potheap.size() == 0) break; // that's all there is
pot = *((gcr->Ix)->pots + (gcr->potheap.front().pot));
// anyway remove top from heap
std::pop_heap(gcr->potheap.begin(), gcr->potheap.end(), hpotcompare);
gcr->potheap.pop_back();
if (pot.LorLeaf == 0) {
// leaf pot - put all the points into the points heap
for (i = 0; i < pot.RorPoints; i++) {
j = pot.points[i];
ct = ((gcr->Ix)->gc + j);
hs.slot = j;
hs.snmd = GeoSNMD(&(gcr->gd), ct);
gcr->slotheap.push_back(hs);
std::push_heap(gcr->slotheap.begin(), gcr->slotheap.end(),
hslotcompare);
}
if (gcr->slotheap.size() != 0) {
slox = gcr->slotheap.front().slot;
gcr->slotsnmd = GeoSNMD(&gcr->gd, (gcr->Ix)->gc + slox);
}
} else {
hp.pot = pot.LorLeaf;
hp.dist = makedist((gcr->Ix)->pots + pot.LorLeaf, &(gcr->gd));
gcr->potheap.push_back(hp);
std::push_heap(gcr->potheap.begin(), gcr->potheap.end(), hpotcompare);
hp.pot = pot.RorPoints;
hp.dist = makedist((gcr->Ix)->pots + pot.RorPoints, &(gcr->gd));
gcr->potheap.push_back(hp);
std::push_heap(gcr->potheap.begin(), gcr->potheap.end(), hpotcompare);
}
gcr->potsnmd = 10.0;
if (gcr->potheap.size() != 0) {
pot = *((gcr->Ix)->pots + (gcr->potheap.front().pot));
gcr->potsnmd = GeoFixtoSNMD(makedist(&pot, &(gcr->gd)));
}
} else {
if (gcr->slotheap.size() == 0) break; // that's all there is
slox = gcr->slotheap.front().slot;
tsnmd = GeoSNMD(&gcr->gd, (gcr->Ix)->gc + slox);
r = GeoResultsGrow(gr);
if (r == -1) {
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gr->snmd);
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gr->slot);
TRI_Free(TRI_UNKNOWN_MEM_ZONE, gr);
return NULL;
}
gr->slot[gr->pointsct] = slox;
gr->snmd[gr->pointsct] = tsnmd;
gr->pointsct++;
gcr->slotsnmd = 5.0;
std::pop_heap(gcr->slotheap.begin(), gcr->slotheap.end(), hslotcompare);
gcr->slotheap.pop_back();
if (gcr->slotheap.size() != 0) {
slox = gcr->slotheap.front().slot;
gcr->slotsnmd = GeoSNMD(&gcr->gd, (gcr->Ix)->gc + slox);
}
}
}
gcts = GeoAnswers(gcr->Ix, gr);
return gcts;
}
void GeoIndex_CursorFree(GeoCursor* gc) {
if (gc == nullptr) {
return;
}
GeoCr* cr = reinterpret_cast<GeoCr*>(gc);
delete cr;
}
/* =================================================== */
/* The remaining routines are usually */
/* only compiled in for debugging purposes. They allow*/
/* the dumping of the index (to a specified file) and */
/* a self-check to see whether the index itself seems */
/* to be correct. */
/* =================================================== */
#ifdef TRI_GEO_DEBUG
void RecursivePotDump(GeoIx* gix, FILE* f, int pot) {
int i;
GeoPot* gp;
gp = gix->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; i < GeoIndexFIXEDPOINTS; i++) fprintf(f, " %x", gp->maxdist[i]);
fprintf(f, "\n");
if (gp->LorLeaf == 0) {
fprintf(f, "Leaf pot containing %d points . . .\n", gp->RorPoints);
for (i = 0; i < gp->RorPoints; i++) {
fprintf(f, "Child %d Point %d ", i, gp->points[i]);
GeoCoordinate* gc = gix->gc + gp->points[i];
fprintf(f, "Lat. %9.4f, Long. %9.4f", gc->latitude, gc->longitude);
#if TRI_GEO_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(gix, f, gp->LorLeaf);
fprintf(f, "\nPot %d - Right Child of pot %d\n", gp->RorPoints, pot);
RecursivePotDump(gix, f, gp->RorPoints);
}
}
void GeoIndex_INDEXDUMP(GeoIdx* gi, FILE* f) {
GeoIx* gix;
gix = (GeoIx*)gi;
fprintf(f, "Dump of entire index. %d pots and %d slots allocated\n",
gix->potct, gix->slotct);
RecursivePotDump(gix, f, 1);
}
int RecursivePotValidate(GeoIx* gix, int pot, int* usage) {
int i, j;
GeoPot* gp;
GeoDetailedPoint gd;
GeoFix maxdist[GeoIndexFIXEDPOINTS];
GeoPot *gpa, *gpb;
gp = gix->pots + pot;
usage[0]++;
if (gp->LorLeaf == 0) {
if ((pot != 1) && (2 * gp->RorPoints < GeoIndexPOTSIZE)) return 1;
for (j = 0; j < GeoIndexFIXEDPOINTS; j++) maxdist[j] = 0;
if (gp->level != 1) return 10;
for (i = 0; i < gp->RorPoints; i++) {
GeoMkDetail(gix, &gd, gix->gc + gp->points[i]);
for (j = 0; j < GeoIndexFIXEDPOINTS; j++)
if (maxdist[j] < gd.fixdist[j]) maxdist[j] = gd.fixdist[j];
if (gd.gs < gp->start) return 8;
if (gd.gs > gp->end) return 9;
}
for (j = 0; j < GeoIndexFIXEDPOINTS; j++)
if (maxdist[j] != gp->maxdist[j]) return 7;
usage[1] += gp->RorPoints;
return 0;
} else {
int pota = gp->LorLeaf;
int potb = gp->RorPoints;
gpa = gix->pots + pota;
gpb = gix->pots + potb;
int lev = gp->level;
int leva = gpa->level;
int levb = gpb->level;
if (leva >= lev) return 2;
if (levb >= lev) return 3;
i = (lev - leva) * (lev - levb);
if (i > 2) return 4;
if (gp->middle != gpa->end) return 5;
if (gp->middle != gpb->start) return 6;
if (gp->start != gpa->start) return 11;
if (gp->end != gpb->end) return 12;
for (j = 0; j < GeoIndexFIXEDPOINTS; j++) maxdist[j] = gpa->maxdist[j];
for (j = 0; j < GeoIndexFIXEDPOINTS; j++)
if (maxdist[j] < gpb->maxdist[j]) maxdist[j] = gpb->maxdist[j];
for (j = 0; j < GeoIndexFIXEDPOINTS; j++)
if (maxdist[j] != gp->maxdist[j]) return 13;
i = RecursivePotValidate(gix, gp->LorLeaf, usage);
if (i != 0) return i;
i = RecursivePotValidate(gix, gp->RorPoints, usage);
if (i != 0) return i;
return 0;
}
}
int GeoIndex_INDEXVALID(GeoIdx* gi) {
int usage[2]; // pots and slots
int j, pot, slot;
GeoIx* gix;
GeoPot* gp;
gix = (GeoIx*)gi;
usage[0] = 0;
usage[1] = 0;
j = RecursivePotValidate(gix, 1, usage);
if (j != 0) return j;
pot = 0;
gp = gix->pots + pot;
pot = gp->LorLeaf;
usage[0]++;
while (pot != 0) {
gp = gix->pots + pot;
pot = gp->LorLeaf;
usage[0]++;
}
if (usage[0] != gix->potct) return 14;
gp = gix->pots + 1;
if (gp->start != 0) return 15;
if (gp->end != 0x1FFFFFFFFFFFFFll) return 16;
slot = 0;
usage[1]++;
slot = (int)((gix->gc[slot]).latitude);
while (slot != 0) {
usage[1]++;
slot = (int)((gix->gc[slot]).latitude);
}
if (usage[1] != gix->slotct) return 17;
return 0;
}
#endif
size_t GeoIndex_MemoryUsage(void* theIndex) {
GeoIx* geoIndex = (GeoIx*)theIndex;
if (geoIndex != NULL) {
return geoIndex->_memoryUsed;
}
return 0;
}
/* end of GeoIndex.c */
Reference in New Issue View Git Blame Copy Permalink