summaryrefslogtreecommitdiffstats
path: root/ext/rtree/geopoly.c
diff options
context:
space:
mode:
Diffstat (limited to 'ext/rtree/geopoly.c')
-rw-r--r--ext/rtree/geopoly.c1814
1 files changed, 1814 insertions, 0 deletions
diff --git a/ext/rtree/geopoly.c b/ext/rtree/geopoly.c
new file mode 100644
index 0000000..7b41e79
--- /dev/null
+++ b/ext/rtree/geopoly.c
@@ -0,0 +1,1814 @@
+/*
+** 2018-05-25
+**
+** The author disclaims copyright to this source code. In place of
+** a legal notice, here is a blessing:
+**
+** May you do good and not evil.
+** May you find forgiveness for yourself and forgive others.
+** May you share freely, never taking more than you give.
+**
+******************************************************************************
+**
+** This file implements an alternative R-Tree virtual table that
+** uses polygons to express the boundaries of 2-dimensional objects.
+**
+** This file is #include-ed onto the end of "rtree.c" so that it has
+** access to all of the R-Tree internals.
+*/
+#include <stdlib.h>
+
+/* Enable -DGEOPOLY_ENABLE_DEBUG for debugging facilities */
+#ifdef GEOPOLY_ENABLE_DEBUG
+ static int geo_debug = 0;
+# define GEODEBUG(X) if(geo_debug)printf X
+#else
+# define GEODEBUG(X)
+#endif
+
+/* Character class routines */
+#ifdef sqlite3Isdigit
+ /* Use the SQLite core versions if this routine is part of the
+ ** SQLite amalgamation */
+# define safe_isdigit(x) sqlite3Isdigit(x)
+# define safe_isalnum(x) sqlite3Isalnum(x)
+# define safe_isxdigit(x) sqlite3Isxdigit(x)
+#else
+ /* Use the standard library for separate compilation */
+#include <ctype.h> /* amalgamator: keep */
+# define safe_isdigit(x) isdigit((unsigned char)(x))
+# define safe_isalnum(x) isalnum((unsigned char)(x))
+# define safe_isxdigit(x) isxdigit((unsigned char)(x))
+#endif
+
+#ifndef JSON_NULL /* The following stuff repeats things found in json1 */
+/*
+** Growing our own isspace() routine this way is twice as fast as
+** the library isspace() function.
+*/
+static const char geopolyIsSpace[] = {
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+};
+#define fast_isspace(x) (geopolyIsSpace[(unsigned char)x])
+#endif /* JSON NULL - back to original code */
+
+/* Compiler and version */
+#ifndef GCC_VERSION
+#if defined(__GNUC__) && !defined(SQLITE_DISABLE_INTRINSIC)
+# define GCC_VERSION (__GNUC__*1000000+__GNUC_MINOR__*1000+__GNUC_PATCHLEVEL__)
+#else
+# define GCC_VERSION 0
+#endif
+#endif
+#ifndef MSVC_VERSION
+#if defined(_MSC_VER) && !defined(SQLITE_DISABLE_INTRINSIC)
+# define MSVC_VERSION _MSC_VER
+#else
+# define MSVC_VERSION 0
+#endif
+#endif
+
+/* Datatype for coordinates
+*/
+typedef float GeoCoord;
+
+/*
+** Internal representation of a polygon.
+**
+** The polygon consists of a sequence of vertexes. There is a line
+** segment between each pair of vertexes, and one final segment from
+** the last vertex back to the first. (This differs from the GeoJSON
+** standard in which the final vertex is a repeat of the first.)
+**
+** The polygon follows the right-hand rule. The area to the right of
+** each segment is "outside" and the area to the left is "inside".
+**
+** The on-disk representation consists of a 4-byte header followed by
+** the values. The 4-byte header is:
+**
+** encoding (1 byte) 0=big-endian, 1=little-endian
+** nvertex (3 bytes) Number of vertexes as a big-endian integer
+**
+** Enough space is allocated for 4 coordinates, to work around over-zealous
+** warnings coming from some compiler (notably, clang). In reality, the size
+** of each GeoPoly memory allocate is adjusted as necessary so that the
+** GeoPoly.a[] array at the end is the appropriate size.
+*/
+typedef struct GeoPoly GeoPoly;
+struct GeoPoly {
+ int nVertex; /* Number of vertexes */
+ unsigned char hdr[4]; /* Header for on-disk representation */
+ GeoCoord a[8]; /* 2*nVertex values. X (longitude) first, then Y */
+};
+
+/* The size of a memory allocation needed for a GeoPoly object sufficient
+** to hold N coordinate pairs.
+*/
+#define GEOPOLY_SZ(N) (sizeof(GeoPoly) + sizeof(GeoCoord)*2*((N)-4))
+
+/* Macros to access coordinates of a GeoPoly.
+** We have to use these macros, rather than just say p->a[i] in order
+** to silence (incorrect) UBSAN warnings if the array index is too large.
+*/
+#define GeoX(P,I) (((GeoCoord*)(P)->a)[(I)*2])
+#define GeoY(P,I) (((GeoCoord*)(P)->a)[(I)*2+1])
+
+
+/*
+** State of a parse of a GeoJSON input.
+*/
+typedef struct GeoParse GeoParse;
+struct GeoParse {
+ const unsigned char *z; /* Unparsed input */
+ int nVertex; /* Number of vertexes in a[] */
+ int nAlloc; /* Space allocated to a[] */
+ int nErr; /* Number of errors encountered */
+ GeoCoord *a; /* Array of vertexes. From sqlite3_malloc64() */
+};
+
+/* Do a 4-byte byte swap */
+static void geopolySwab32(unsigned char *a){
+ unsigned char t = a[0];
+ a[0] = a[3];
+ a[3] = t;
+ t = a[1];
+ a[1] = a[2];
+ a[2] = t;
+}
+
+/* Skip whitespace. Return the next non-whitespace character. */
+static char geopolySkipSpace(GeoParse *p){
+ while( fast_isspace(p->z[0]) ) p->z++;
+ return p->z[0];
+}
+
+/* Parse out a number. Write the value into *pVal if pVal!=0.
+** return non-zero on success and zero if the next token is not a number.
+*/
+static int geopolyParseNumber(GeoParse *p, GeoCoord *pVal){
+ char c = geopolySkipSpace(p);
+ const unsigned char *z = p->z;
+ int j = 0;
+ int seenDP = 0;
+ int seenE = 0;
+ if( c=='-' ){
+ j = 1;
+ c = z[j];
+ }
+ if( c=='0' && z[j+1]>='0' && z[j+1]<='9' ) return 0;
+ for(;; j++){
+ c = z[j];
+ if( safe_isdigit(c) ) continue;
+ if( c=='.' ){
+ if( z[j-1]=='-' ) return 0;
+ if( seenDP ) return 0;
+ seenDP = 1;
+ continue;
+ }
+ if( c=='e' || c=='E' ){
+ if( z[j-1]<'0' ) return 0;
+ if( seenE ) return -1;
+ seenDP = seenE = 1;
+ c = z[j+1];
+ if( c=='+' || c=='-' ){
+ j++;
+ c = z[j+1];
+ }
+ if( c<'0' || c>'9' ) return 0;
+ continue;
+ }
+ break;
+ }
+ if( z[j-1]<'0' ) return 0;
+ if( pVal ){
+#ifdef SQLITE_AMALGAMATION
+ /* The sqlite3AtoF() routine is much much faster than atof(), if it
+ ** is available */
+ double r;
+ (void)sqlite3AtoF((const char*)p->z, &r, j, SQLITE_UTF8);
+ *pVal = r;
+#else
+ *pVal = (GeoCoord)atof((const char*)p->z);
+#endif
+ }
+ p->z += j;
+ return 1;
+}
+
+/*
+** If the input is a well-formed JSON array of coordinates with at least
+** four coordinates and where each coordinate is itself a two-value array,
+** then convert the JSON into a GeoPoly object and return a pointer to
+** that object.
+**
+** If any error occurs, return NULL.
+*/
+static GeoPoly *geopolyParseJson(const unsigned char *z, int *pRc){
+ GeoParse s;
+ int rc = SQLITE_OK;
+ memset(&s, 0, sizeof(s));
+ s.z = z;
+ if( geopolySkipSpace(&s)=='[' ){
+ s.z++;
+ while( geopolySkipSpace(&s)=='[' ){
+ int ii = 0;
+ char c;
+ s.z++;
+ if( s.nVertex>=s.nAlloc ){
+ GeoCoord *aNew;
+ s.nAlloc = s.nAlloc*2 + 16;
+ aNew = sqlite3_realloc64(s.a, s.nAlloc*sizeof(GeoCoord)*2 );
+ if( aNew==0 ){
+ rc = SQLITE_NOMEM;
+ s.nErr++;
+ break;
+ }
+ s.a = aNew;
+ }
+ while( geopolyParseNumber(&s, ii<=1 ? &s.a[s.nVertex*2+ii] : 0) ){
+ ii++;
+ if( ii==2 ) s.nVertex++;
+ c = geopolySkipSpace(&s);
+ s.z++;
+ if( c==',' ) continue;
+ if( c==']' && ii>=2 ) break;
+ s.nErr++;
+ rc = SQLITE_ERROR;
+ goto parse_json_err;
+ }
+ if( geopolySkipSpace(&s)==',' ){
+ s.z++;
+ continue;
+ }
+ break;
+ }
+ if( geopolySkipSpace(&s)==']'
+ && s.nVertex>=4
+ && s.a[0]==s.a[s.nVertex*2-2]
+ && s.a[1]==s.a[s.nVertex*2-1]
+ && (s.z++, geopolySkipSpace(&s)==0)
+ ){
+ GeoPoly *pOut;
+ int x = 1;
+ s.nVertex--; /* Remove the redundant vertex at the end */
+ pOut = sqlite3_malloc64( GEOPOLY_SZ((sqlite3_int64)s.nVertex) );
+ x = 1;
+ if( pOut==0 ) goto parse_json_err;
+ pOut->nVertex = s.nVertex;
+ memcpy(pOut->a, s.a, s.nVertex*2*sizeof(GeoCoord));
+ pOut->hdr[0] = *(unsigned char*)&x;
+ pOut->hdr[1] = (s.nVertex>>16)&0xff;
+ pOut->hdr[2] = (s.nVertex>>8)&0xff;
+ pOut->hdr[3] = s.nVertex&0xff;
+ sqlite3_free(s.a);
+ if( pRc ) *pRc = SQLITE_OK;
+ return pOut;
+ }else{
+ s.nErr++;
+ rc = SQLITE_ERROR;
+ }
+ }
+parse_json_err:
+ if( pRc ) *pRc = rc;
+ sqlite3_free(s.a);
+ return 0;
+}
+
+/*
+** Given a function parameter, try to interpret it as a polygon, either
+** in the binary format or JSON text. Compute a GeoPoly object and
+** return a pointer to that object. Or if the input is not a well-formed
+** polygon, put an error message in sqlite3_context and return NULL.
+*/
+static GeoPoly *geopolyFuncParam(
+ sqlite3_context *pCtx, /* Context for error messages */
+ sqlite3_value *pVal, /* The value to decode */
+ int *pRc /* Write error here */
+){
+ GeoPoly *p = 0;
+ int nByte;
+ testcase( pCtx==0 );
+ if( sqlite3_value_type(pVal)==SQLITE_BLOB
+ && (nByte = sqlite3_value_bytes(pVal))>=(4+6*sizeof(GeoCoord))
+ ){
+ const unsigned char *a = sqlite3_value_blob(pVal);
+ int nVertex;
+ if( a==0 ){
+ if( pCtx ) sqlite3_result_error_nomem(pCtx);
+ return 0;
+ }
+ nVertex = (a[1]<<16) + (a[2]<<8) + a[3];
+ if( (a[0]==0 || a[0]==1)
+ && (nVertex*2*sizeof(GeoCoord) + 4)==(unsigned int)nByte
+ ){
+ p = sqlite3_malloc64( sizeof(*p) + (nVertex-1)*2*sizeof(GeoCoord) );
+ if( p==0 ){
+ if( pRc ) *pRc = SQLITE_NOMEM;
+ if( pCtx ) sqlite3_result_error_nomem(pCtx);
+ }else{
+ int x = 1;
+ p->nVertex = nVertex;
+ memcpy(p->hdr, a, nByte);
+ if( a[0] != *(unsigned char*)&x ){
+ int ii;
+ for(ii=0; ii<nVertex; ii++){
+ geopolySwab32((unsigned char*)&GeoX(p,ii));
+ geopolySwab32((unsigned char*)&GeoY(p,ii));
+ }
+ p->hdr[0] ^= 1;
+ }
+ }
+ }
+ if( pRc ) *pRc = SQLITE_OK;
+ return p;
+ }else if( sqlite3_value_type(pVal)==SQLITE_TEXT ){
+ const unsigned char *zJson = sqlite3_value_text(pVal);
+ if( zJson==0 ){
+ if( pRc ) *pRc = SQLITE_NOMEM;
+ return 0;
+ }
+ return geopolyParseJson(zJson, pRc);
+ }else{
+ if( pRc ) *pRc = SQLITE_ERROR;
+ return 0;
+ }
+}
+
+/*
+** Implementation of the geopoly_blob(X) function.
+**
+** If the input is a well-formed Geopoly BLOB or JSON string
+** then return the BLOB representation of the polygon. Otherwise
+** return NULL.
+*/
+static void geopolyBlobFunc(
+ sqlite3_context *context,
+ int argc,
+ sqlite3_value **argv
+){
+ GeoPoly *p = geopolyFuncParam(context, argv[0], 0);
+ if( p ){
+ sqlite3_result_blob(context, p->hdr,
+ 4+8*p->nVertex, SQLITE_TRANSIENT);
+ sqlite3_free(p);
+ }
+}
+
+/*
+** SQL function: geopoly_json(X)
+**
+** Interpret X as a polygon and render it as a JSON array
+** of coordinates. Or, if X is not a valid polygon, return NULL.
+*/
+static void geopolyJsonFunc(
+ sqlite3_context *context,
+ int argc,
+ sqlite3_value **argv
+){
+ GeoPoly *p = geopolyFuncParam(context, argv[0], 0);
+ if( p ){
+ sqlite3 *db = sqlite3_context_db_handle(context);
+ sqlite3_str *x = sqlite3_str_new(db);
+ int i;
+ sqlite3_str_append(x, "[", 1);
+ for(i=0; i<p->nVertex; i++){
+ sqlite3_str_appendf(x, "[%!g,%!g],", GeoX(p,i), GeoY(p,i));
+ }
+ sqlite3_str_appendf(x, "[%!g,%!g]]", GeoX(p,0), GeoY(p,0));
+ sqlite3_result_text(context, sqlite3_str_finish(x), -1, sqlite3_free);
+ sqlite3_free(p);
+ }
+}
+
+/*
+** SQL function: geopoly_svg(X, ....)
+**
+** Interpret X as a polygon and render it as a SVG <polyline>.
+** Additional arguments are added as attributes to the <polyline>.
+*/
+static void geopolySvgFunc(
+ sqlite3_context *context,
+ int argc,
+ sqlite3_value **argv
+){
+ GeoPoly *p;
+ if( argc<1 ) return;
+ p = geopolyFuncParam(context, argv[0], 0);
+ if( p ){
+ sqlite3 *db = sqlite3_context_db_handle(context);
+ sqlite3_str *x = sqlite3_str_new(db);
+ int i;
+ char cSep = '\'';
+ sqlite3_str_appendf(x, "<polyline points=");
+ for(i=0; i<p->nVertex; i++){
+ sqlite3_str_appendf(x, "%c%g,%g", cSep, GeoX(p,i), GeoY(p,i));
+ cSep = ' ';
+ }
+ sqlite3_str_appendf(x, " %g,%g'", GeoX(p,0), GeoY(p,0));
+ for(i=1; i<argc; i++){
+ const char *z = (const char*)sqlite3_value_text(argv[i]);
+ if( z && z[0] ){
+ sqlite3_str_appendf(x, " %s", z);
+ }
+ }
+ sqlite3_str_appendf(x, "></polyline>");
+ sqlite3_result_text(context, sqlite3_str_finish(x), -1, sqlite3_free);
+ sqlite3_free(p);
+ }
+}
+
+/*
+** SQL Function: geopoly_xform(poly, A, B, C, D, E, F)
+**
+** Transform and/or translate a polygon as follows:
+**
+** x1 = A*x0 + B*y0 + E
+** y1 = C*x0 + D*y0 + F
+**
+** For a translation:
+**
+** geopoly_xform(poly, 1, 0, 0, 1, x-offset, y-offset)
+**
+** Rotate by R around the point (0,0):
+**
+** geopoly_xform(poly, cos(R), sin(R), -sin(R), cos(R), 0, 0)
+*/
+static void geopolyXformFunc(
+ sqlite3_context *context,
+ int argc,
+ sqlite3_value **argv
+){
+ GeoPoly *p = geopolyFuncParam(context, argv[0], 0);
+ double A = sqlite3_value_double(argv[1]);
+ double B = sqlite3_value_double(argv[2]);
+ double C = sqlite3_value_double(argv[3]);
+ double D = sqlite3_value_double(argv[4]);
+ double E = sqlite3_value_double(argv[5]);
+ double F = sqlite3_value_double(argv[6]);
+ GeoCoord x1, y1, x0, y0;
+ int ii;
+ if( p ){
+ for(ii=0; ii<p->nVertex; ii++){
+ x0 = GeoX(p,ii);
+ y0 = GeoY(p,ii);
+ x1 = (GeoCoord)(A*x0 + B*y0 + E);
+ y1 = (GeoCoord)(C*x0 + D*y0 + F);
+ GeoX(p,ii) = x1;
+ GeoY(p,ii) = y1;
+ }
+ sqlite3_result_blob(context, p->hdr,
+ 4+8*p->nVertex, SQLITE_TRANSIENT);
+ sqlite3_free(p);
+ }
+}
+
+/*
+** Compute the area enclosed by the polygon.
+**
+** This routine can also be used to detect polygons that rotate in
+** the wrong direction. Polygons are suppose to be counter-clockwise (CCW).
+** This routine returns a negative value for clockwise (CW) polygons.
+*/
+static double geopolyArea(GeoPoly *p){
+ double rArea = 0.0;
+ int ii;
+ for(ii=0; ii<p->nVertex-1; ii++){
+ rArea += (GeoX(p,ii) - GeoX(p,ii+1)) /* (x0 - x1) */
+ * (GeoY(p,ii) + GeoY(p,ii+1)) /* (y0 + y1) */
+ * 0.5;
+ }
+ rArea += (GeoX(p,ii) - GeoX(p,0)) /* (xN - x0) */
+ * (GeoY(p,ii) + GeoY(p,0)) /* (yN + y0) */
+ * 0.5;
+ return rArea;
+}
+
+/*
+** Implementation of the geopoly_area(X) function.
+**
+** If the input is a well-formed Geopoly BLOB then return the area
+** enclosed by the polygon. If the polygon circulates clockwise instead
+** of counterclockwise (as it should) then return the negative of the
+** enclosed area. Otherwise return NULL.
+*/
+static void geopolyAreaFunc(
+ sqlite3_context *context,
+ int argc,
+ sqlite3_value **argv
+){
+ GeoPoly *p = geopolyFuncParam(context, argv[0], 0);
+ if( p ){
+ sqlite3_result_double(context, geopolyArea(p));
+ sqlite3_free(p);
+ }
+}
+
+/*
+** Implementation of the geopoly_ccw(X) function.
+**
+** If the rotation of polygon X is clockwise (incorrect) instead of
+** counter-clockwise (the correct winding order according to RFC7946)
+** then reverse the order of the vertexes in polygon X.
+**
+** In other words, this routine returns a CCW polygon regardless of the
+** winding order of its input.
+**
+** Use this routine to sanitize historical inputs that that sometimes
+** contain polygons that wind in the wrong direction.
+*/
+static void geopolyCcwFunc(
+ sqlite3_context *context,
+ int argc,
+ sqlite3_value **argv
+){
+ GeoPoly *p = geopolyFuncParam(context, argv[0], 0);
+ if( p ){
+ if( geopolyArea(p)<0.0 ){
+ int ii, jj;
+ for(ii=1, jj=p->nVertex-1; ii<jj; ii++, jj--){
+ GeoCoord t = GeoX(p,ii);
+ GeoX(p,ii) = GeoX(p,jj);
+ GeoX(p,jj) = t;
+ t = GeoY(p,ii);
+ GeoY(p,ii) = GeoY(p,jj);
+ GeoY(p,jj) = t;
+ }
+ }
+ sqlite3_result_blob(context, p->hdr,
+ 4+8*p->nVertex, SQLITE_TRANSIENT);
+ sqlite3_free(p);
+ }
+}
+
+#define GEOPOLY_PI 3.1415926535897932385
+
+/* Fast approximation for sine(X) for X between -0.5*pi and 2*pi
+*/
+static double geopolySine(double r){
+ assert( r>=-0.5*GEOPOLY_PI && r<=2.0*GEOPOLY_PI );
+ if( r>=1.5*GEOPOLY_PI ){
+ r -= 2.0*GEOPOLY_PI;
+ }
+ if( r>=0.5*GEOPOLY_PI ){
+ return -geopolySine(r-GEOPOLY_PI);
+ }else{
+ double r2 = r*r;
+ double r3 = r2*r;
+ double r5 = r3*r2;
+ return 0.9996949*r - 0.1656700*r3 + 0.0075134*r5;
+ }
+}
+
+/*
+** Function: geopoly_regular(X,Y,R,N)
+**
+** Construct a simple, convex, regular polygon centered at X, Y
+** with circumradius R and with N sides.
+*/
+static void geopolyRegularFunc(
+ sqlite3_context *context,
+ int argc,
+ sqlite3_value **argv
+){
+ double x = sqlite3_value_double(argv[0]);
+ double y = sqlite3_value_double(argv[1]);
+ double r = sqlite3_value_double(argv[2]);
+ int n = sqlite3_value_int(argv[3]);
+ int i;
+ GeoPoly *p;
+
+ if( n<3 || r<=0.0 ) return;
+ if( n>1000 ) n = 1000;
+ p = sqlite3_malloc64( sizeof(*p) + (n-1)*2*sizeof(GeoCoord) );
+ if( p==0 ){
+ sqlite3_result_error_nomem(context);
+ return;
+ }
+ i = 1;
+ p->hdr[0] = *(unsigned char*)&i;
+ p->hdr[1] = 0;
+ p->hdr[2] = (n>>8)&0xff;
+ p->hdr[3] = n&0xff;
+ for(i=0; i<n; i++){
+ double rAngle = 2.0*GEOPOLY_PI*i/n;
+ GeoX(p,i) = x - r*geopolySine(rAngle-0.5*GEOPOLY_PI);
+ GeoY(p,i) = y + r*geopolySine(rAngle);
+ }
+ sqlite3_result_blob(context, p->hdr, 4+8*n, SQLITE_TRANSIENT);
+ sqlite3_free(p);
+}
+
+/*
+** If pPoly is a polygon, compute its bounding box. Then:
+**
+** (1) if aCoord!=0 store the bounding box in aCoord, returning NULL
+** (2) otherwise, compute a GeoPoly for the bounding box and return the
+** new GeoPoly
+**
+** If pPoly is NULL but aCoord is not NULL, then compute a new GeoPoly from
+** the bounding box in aCoord and return a pointer to that GeoPoly.
+*/
+static GeoPoly *geopolyBBox(
+ sqlite3_context *context, /* For recording the error */
+ sqlite3_value *pPoly, /* The polygon */
+ RtreeCoord *aCoord, /* Results here */
+ int *pRc /* Error code here */
+){
+ GeoPoly *pOut = 0;
+ GeoPoly *p;
+ float mnX, mxX, mnY, mxY;
+ if( pPoly==0 && aCoord!=0 ){
+ p = 0;
+ mnX = aCoord[0].f;
+ mxX = aCoord[1].f;
+ mnY = aCoord[2].f;
+ mxY = aCoord[3].f;
+ goto geopolyBboxFill;
+ }else{
+ p = geopolyFuncParam(context, pPoly, pRc);
+ }
+ if( p ){
+ int ii;
+ mnX = mxX = GeoX(p,0);
+ mnY = mxY = GeoY(p,0);
+ for(ii=1; ii<p->nVertex; ii++){
+ double r = GeoX(p,ii);
+ if( r<mnX ) mnX = (float)r;
+ else if( r>mxX ) mxX = (float)r;
+ r = GeoY(p,ii);
+ if( r<mnY ) mnY = (float)r;
+ else if( r>mxY ) mxY = (float)r;
+ }
+ if( pRc ) *pRc = SQLITE_OK;
+ if( aCoord==0 ){
+ geopolyBboxFill:
+ pOut = sqlite3_realloc64(p, GEOPOLY_SZ(4));
+ if( pOut==0 ){
+ sqlite3_free(p);
+ if( context ) sqlite3_result_error_nomem(context);
+ if( pRc ) *pRc = SQLITE_NOMEM;
+ return 0;
+ }
+ pOut->nVertex = 4;
+ ii = 1;
+ pOut->hdr[0] = *(unsigned char*)&ii;
+ pOut->hdr[1] = 0;
+ pOut->hdr[2] = 0;
+ pOut->hdr[3] = 4;
+ GeoX(pOut,0) = mnX;
+ GeoY(pOut,0) = mnY;
+ GeoX(pOut,1) = mxX;
+ GeoY(pOut,1) = mnY;
+ GeoX(pOut,2) = mxX;
+ GeoY(pOut,2) = mxY;
+ GeoX(pOut,3) = mnX;
+ GeoY(pOut,3) = mxY;
+ }else{
+ sqlite3_free(p);
+ aCoord[0].f = mnX;
+ aCoord[1].f = mxX;
+ aCoord[2].f = mnY;
+ aCoord[3].f = mxY;
+ }
+ }else if( aCoord ){
+ memset(aCoord, 0, sizeof(RtreeCoord)*4);
+ }
+ return pOut;
+}
+
+/*
+** Implementation of the geopoly_bbox(X) SQL function.
+*/
+static void geopolyBBoxFunc(
+ sqlite3_context *context,
+ int argc,
+ sqlite3_value **argv
+){
+ GeoPoly *p = geopolyBBox(context, argv[0], 0, 0);
+ if( p ){
+ sqlite3_result_blob(context, p->hdr,
+ 4+8*p->nVertex, SQLITE_TRANSIENT);
+ sqlite3_free(p);
+ }
+}
+
+/*
+** State vector for the geopoly_group_bbox() aggregate function.
+*/
+typedef struct GeoBBox GeoBBox;
+struct GeoBBox {
+ int isInit;
+ RtreeCoord a[4];
+};
+
+
+/*
+** Implementation of the geopoly_group_bbox(X) aggregate SQL function.
+*/
+static void geopolyBBoxStep(
+ sqlite3_context *context,
+ int argc,
+ sqlite3_value **argv
+){
+ RtreeCoord a[4];
+ int rc = SQLITE_OK;
+ (void)geopolyBBox(context, argv[0], a, &rc);
+ if( rc==SQLITE_OK ){
+ GeoBBox *pBBox;
+ pBBox = (GeoBBox*)sqlite3_aggregate_context(context, sizeof(*pBBox));
+ if( pBBox==0 ) return;
+ if( pBBox->isInit==0 ){
+ pBBox->isInit = 1;
+ memcpy(pBBox->a, a, sizeof(RtreeCoord)*4);
+ }else{
+ if( a[0].f < pBBox->a[0].f ) pBBox->a[0] = a[0];
+ if( a[1].f > pBBox->a[1].f ) pBBox->a[1] = a[1];
+ if( a[2].f < pBBox->a[2].f ) pBBox->a[2] = a[2];
+ if( a[3].f > pBBox->a[3].f ) pBBox->a[3] = a[3];
+ }
+ }
+}
+static void geopolyBBoxFinal(
+ sqlite3_context *context
+){
+ GeoPoly *p;
+ GeoBBox *pBBox;
+ pBBox = (GeoBBox*)sqlite3_aggregate_context(context, 0);
+ if( pBBox==0 ) return;
+ p = geopolyBBox(context, 0, pBBox->a, 0);
+ if( p ){
+ sqlite3_result_blob(context, p->hdr,
+ 4+8*p->nVertex, SQLITE_TRANSIENT);
+ sqlite3_free(p);
+ }
+}
+
+
+/*
+** Determine if point (x0,y0) is beneath line segment (x1,y1)->(x2,y2).
+** Returns:
+**
+** +2 x0,y0 is on the line segement
+**
+** +1 x0,y0 is beneath line segment
+**
+** 0 x0,y0 is not on or beneath the line segment or the line segment
+** is vertical and x0,y0 is not on the line segment
+**
+** The left-most coordinate min(x1,x2) is not considered to be part of
+** the line segment for the purposes of this analysis.
+*/
+static int pointBeneathLine(
+ double x0, double y0,
+ double x1, double y1,
+ double x2, double y2
+){
+ double y;
+ if( x0==x1 && y0==y1 ) return 2;
+ if( x1<x2 ){
+ if( x0<=x1 || x0>x2 ) return 0;
+ }else if( x1>x2 ){
+ if( x0<=x2 || x0>x1 ) return 0;
+ }else{
+ /* Vertical line segment */
+ if( x0!=x1 ) return 0;
+ if( y0<y1 && y0<y2 ) return 0;
+ if( y0>y1 && y0>y2 ) return 0;
+ return 2;
+ }
+ y = y1 + (y2-y1)*(x0-x1)/(x2-x1);
+ if( y0==y ) return 2;
+ if( y0<y ) return 1;
+ return 0;
+}
+
+/*
+** SQL function: geopoly_contains_point(P,X,Y)
+**
+** Return +2 if point X,Y is within polygon P.
+** Return +1 if point X,Y is on the polygon boundary.
+** Return 0 if point X,Y is outside the polygon
+*/
+static void geopolyContainsPointFunc(
+ sqlite3_context *context,
+ int argc,
+ sqlite3_value **argv
+){
+ GeoPoly *p1 = geopolyFuncParam(context, argv[0], 0);
+ double x0 = sqlite3_value_double(argv[1]);
+ double y0 = sqlite3_value_double(argv[2]);
+ int v = 0;
+ int cnt = 0;
+ int ii;
+ if( p1==0 ) return;
+ for(ii=0; ii<p1->nVertex-1; ii++){
+ v = pointBeneathLine(x0,y0,GeoX(p1,ii), GeoY(p1,ii),
+ GeoX(p1,ii+1),GeoY(p1,ii+1));
+ if( v==2 ) break;
+ cnt += v;
+ }
+ if( v!=2 ){
+ v = pointBeneathLine(x0,y0,GeoX(p1,ii), GeoY(p1,ii),
+ GeoX(p1,0), GeoY(p1,0));
+ }
+ if( v==2 ){
+ sqlite3_result_int(context, 1);
+ }else if( ((v+cnt)&1)==0 ){
+ sqlite3_result_int(context, 0);
+ }else{
+ sqlite3_result_int(context, 2);
+ }
+ sqlite3_free(p1);
+}
+
+/* Forward declaration */
+static int geopolyOverlap(GeoPoly *p1, GeoPoly *p2);
+
+/*
+** SQL function: geopoly_within(P1,P2)
+**
+** Return +2 if P1 and P2 are the same polygon
+** Return +1 if P2 is contained within P1
+** Return 0 if any part of P2 is on the outside of P1
+**
+*/
+static void geopolyWithinFunc(
+ sqlite3_context *context,
+ int argc,
+ sqlite3_value **argv
+){
+ GeoPoly *p1 = geopolyFuncParam(context, argv[0], 0);
+ GeoPoly *p2 = geopolyFuncParam(context, argv[1], 0);
+ if( p1 && p2 ){
+ int x = geopolyOverlap(p1, p2);
+ if( x<0 ){
+ sqlite3_result_error_nomem(context);
+ }else{
+ sqlite3_result_int(context, x==2 ? 1 : x==4 ? 2 : 0);
+ }
+ }
+ sqlite3_free(p1);
+ sqlite3_free(p2);
+}
+
+/* Objects used by the overlap algorihm. */
+typedef struct GeoEvent GeoEvent;
+typedef struct GeoSegment GeoSegment;
+typedef struct GeoOverlap GeoOverlap;
+struct GeoEvent {
+ double x; /* X coordinate at which event occurs */
+ int eType; /* 0 for ADD, 1 for REMOVE */
+ GeoSegment *pSeg; /* The segment to be added or removed */
+ GeoEvent *pNext; /* Next event in the sorted list */
+};
+struct GeoSegment {
+ double C, B; /* y = C*x + B */
+ double y; /* Current y value */
+ float y0; /* Initial y value */
+ unsigned char side; /* 1 for p1, 2 for p2 */
+ unsigned int idx; /* Which segment within the side */
+ GeoSegment *pNext; /* Next segment in a list sorted by y */
+};
+struct GeoOverlap {
+ GeoEvent *aEvent; /* Array of all events */
+ GeoSegment *aSegment; /* Array of all segments */
+ int nEvent; /* Number of events */
+ int nSegment; /* Number of segments */
+};
+
+/*
+** Add a single segment and its associated events.
+*/
+static void geopolyAddOneSegment(
+ GeoOverlap *p,
+ GeoCoord x0,
+ GeoCoord y0,
+ GeoCoord x1,
+ GeoCoord y1,
+ unsigned char side,
+ unsigned int idx
+){
+ GeoSegment *pSeg;
+ GeoEvent *pEvent;
+ if( x0==x1 ) return; /* Ignore vertical segments */
+ if( x0>x1 ){
+ GeoCoord t = x0;
+ x0 = x1;
+ x1 = t;
+ t = y0;
+ y0 = y1;
+ y1 = t;
+ }
+ pSeg = p->aSegment + p->nSegment;
+ p->nSegment++;
+ pSeg->C = (y1-y0)/(x1-x0);
+ pSeg->B = y1 - x1*pSeg->C;
+ pSeg->y0 = y0;
+ pSeg->side = side;
+ pSeg->idx = idx;
+ pEvent = p->aEvent + p->nEvent;
+ p->nEvent++;
+ pEvent->x = x0;
+ pEvent->eType = 0;
+ pEvent->pSeg = pSeg;
+ pEvent = p->aEvent + p->nEvent;
+ p->nEvent++;
+ pEvent->x = x1;
+ pEvent->eType = 1;
+ pEvent->pSeg = pSeg;
+}
+
+
+
+/*
+** Insert all segments and events for polygon pPoly.
+*/
+static void geopolyAddSegments(
+ GeoOverlap *p, /* Add segments to this Overlap object */
+ GeoPoly *pPoly, /* Take all segments from this polygon */
+ unsigned char side /* The side of pPoly */
+){
+ unsigned int i;
+ GeoCoord *x;
+ for(i=0; i<(unsigned)pPoly->nVertex-1; i++){
+ x = &GeoX(pPoly,i);
+ geopolyAddOneSegment(p, x[0], x[1], x[2], x[3], side, i);
+ }
+ x = &GeoX(pPoly,i);
+ geopolyAddOneSegment(p, x[0], x[1], pPoly->a[0], pPoly->a[1], side, i);
+}
+
+/*
+** Merge two lists of sorted events by X coordinate
+*/
+static GeoEvent *geopolyEventMerge(GeoEvent *pLeft, GeoEvent *pRight){
+ GeoEvent head, *pLast;
+ head.pNext = 0;
+ pLast = &head;
+ while( pRight && pLeft ){
+ if( pRight->x <= pLeft->x ){
+ pLast->pNext = pRight;
+ pLast = pRight;
+ pRight = pRight->pNext;
+ }else{
+ pLast->pNext = pLeft;
+ pLast = pLeft;
+ pLeft = pLeft->pNext;
+ }
+ }
+ pLast->pNext = pRight ? pRight : pLeft;
+ return head.pNext;
+}
+
+/*
+** Sort an array of nEvent event objects into a list.
+*/
+static GeoEvent *geopolySortEventsByX(GeoEvent *aEvent, int nEvent){
+ int mx = 0;
+ int i, j;
+ GeoEvent *p;
+ GeoEvent *a[50];
+ for(i=0; i<nEvent; i++){
+ p = &aEvent[i];
+ p->pNext = 0;
+ for(j=0; j<mx && a[j]; j++){
+ p = geopolyEventMerge(a[j], p);
+ a[j] = 0;
+ }
+ a[j] = p;
+ if( j>=mx ) mx = j+1;
+ }
+ p = 0;
+ for(i=0; i<mx; i++){
+ p = geopolyEventMerge(a[i], p);
+ }
+ return p;
+}
+
+/*
+** Merge two lists of sorted segments by Y, and then by C.
+*/
+static GeoSegment *geopolySegmentMerge(GeoSegment *pLeft, GeoSegment *pRight){
+ GeoSegment head, *pLast;
+ head.pNext = 0;
+ pLast = &head;
+ while( pRight && pLeft ){
+ double r = pRight->y - pLeft->y;
+ if( r==0.0 ) r = pRight->C - pLeft->C;
+ if( r<0.0 ){
+ pLast->pNext = pRight;
+ pLast = pRight;
+ pRight = pRight->pNext;
+ }else{
+ pLast->pNext = pLeft;
+ pLast = pLeft;
+ pLeft = pLeft->pNext;
+ }
+ }
+ pLast->pNext = pRight ? pRight : pLeft;
+ return head.pNext;
+}
+
+/*
+** Sort a list of GeoSegments in order of increasing Y and in the event of
+** a tie, increasing C (slope).
+*/
+static GeoSegment *geopolySortSegmentsByYAndC(GeoSegment *pList){
+ int mx = 0;
+ int i;
+ GeoSegment *p;
+ GeoSegment *a[50];
+ while( pList ){
+ p = pList;
+ pList = pList->pNext;
+ p->pNext = 0;
+ for(i=0; i<mx && a[i]; i++){
+ p = geopolySegmentMerge(a[i], p);
+ a[i] = 0;
+ }
+ a[i] = p;
+ if( i>=mx ) mx = i+1;
+ }
+ p = 0;
+ for(i=0; i<mx; i++){
+ p = geopolySegmentMerge(a[i], p);
+ }
+ return p;
+}
+
+/*
+** Determine the overlap between two polygons
+*/
+static int geopolyOverlap(GeoPoly *p1, GeoPoly *p2){
+ sqlite3_int64 nVertex = p1->nVertex + p2->nVertex + 2;
+ GeoOverlap *p;
+ sqlite3_int64 nByte;
+ GeoEvent *pThisEvent;
+ double rX;
+ int rc = 0;
+ int needSort = 0;
+ GeoSegment *pActive = 0;
+ GeoSegment *pSeg;
+ unsigned char aOverlap[4];
+
+ nByte = sizeof(GeoEvent)*nVertex*2
+ + sizeof(GeoSegment)*nVertex
+ + sizeof(GeoOverlap);
+ p = sqlite3_malloc64( nByte );
+ if( p==0 ) return -1;
+ p->aEvent = (GeoEvent*)&p[1];
+ p->aSegment = (GeoSegment*)&p->aEvent[nVertex*2];
+ p->nEvent = p->nSegment = 0;
+ geopolyAddSegments(p, p1, 1);
+ geopolyAddSegments(p, p2, 2);
+ pThisEvent = geopolySortEventsByX(p->aEvent, p->nEvent);
+ rX = pThisEvent && pThisEvent->x==0.0 ? -1.0 : 0.0;
+ memset(aOverlap, 0, sizeof(aOverlap));
+ while( pThisEvent ){
+ if( pThisEvent->x!=rX ){
+ GeoSegment *pPrev = 0;
+ int iMask = 0;
+ GEODEBUG(("Distinct X: %g\n", pThisEvent->x));
+ rX = pThisEvent->x;
+ if( needSort ){
+ GEODEBUG(("SORT\n"));
+ pActive = geopolySortSegmentsByYAndC(pActive);
+ needSort = 0;
+ }
+ for(pSeg=pActive; pSeg; pSeg=pSeg->pNext){
+ if( pPrev ){
+ if( pPrev->y!=pSeg->y ){
+ GEODEBUG(("MASK: %d\n", iMask));
+ aOverlap[iMask] = 1;
+ }
+ }
+ iMask ^= pSeg->side;
+ pPrev = pSeg;
+ }
+ pPrev = 0;
+ for(pSeg=pActive; pSeg; pSeg=pSeg->pNext){
+ double y = pSeg->C*rX + pSeg->B;
+ GEODEBUG(("Segment %d.%d %g->%g\n", pSeg->side, pSeg->idx, pSeg->y, y));
+ pSeg->y = y;
+ if( pPrev ){
+ if( pPrev->y>pSeg->y && pPrev->side!=pSeg->side ){
+ rc = 1;
+ GEODEBUG(("Crossing: %d.%d and %d.%d\n",
+ pPrev->side, pPrev->idx,
+ pSeg->side, pSeg->idx));
+ goto geopolyOverlapDone;
+ }else if( pPrev->y!=pSeg->y ){
+ GEODEBUG(("MASK: %d\n", iMask));
+ aOverlap[iMask] = 1;
+ }
+ }
+ iMask ^= pSeg->side;
+ pPrev = pSeg;
+ }
+ }
+ GEODEBUG(("%s %d.%d C=%g B=%g\n",
+ pThisEvent->eType ? "RM " : "ADD",
+ pThisEvent->pSeg->side, pThisEvent->pSeg->idx,
+ pThisEvent->pSeg->C,
+ pThisEvent->pSeg->B));
+ if( pThisEvent->eType==0 ){
+ /* Add a segment */
+ pSeg = pThisEvent->pSeg;
+ pSeg->y = pSeg->y0;
+ pSeg->pNext = pActive;
+ pActive = pSeg;
+ needSort = 1;
+ }else{
+ /* Remove a segment */
+ if( pActive==pThisEvent->pSeg ){
+ pActive = ALWAYS(pActive) ? pActive->pNext : 0;
+ }else{
+ for(pSeg=pActive; pSeg; pSeg=pSeg->pNext){
+ if( pSeg->pNext==pThisEvent->pSeg ){
+ pSeg->pNext = ALWAYS(pSeg->pNext) ? pSeg->pNext->pNext : 0;
+ break;
+ }
+ }
+ }
+ }
+ pThisEvent = pThisEvent->pNext;
+ }
+ if( aOverlap[3]==0 ){
+ rc = 0;
+ }else if( aOverlap[1]!=0 && aOverlap[2]==0 ){
+ rc = 3;
+ }else if( aOverlap[1]==0 && aOverlap[2]!=0 ){
+ rc = 2;
+ }else if( aOverlap[1]==0 && aOverlap[2]==0 ){
+ rc = 4;
+ }else{
+ rc = 1;
+ }
+
+geopolyOverlapDone:
+ sqlite3_free(p);
+ return rc;
+}
+
+/*
+** SQL function: geopoly_overlap(P1,P2)
+**
+** Determine whether or not P1 and P2 overlap. Return value:
+**
+** 0 The two polygons are disjoint
+** 1 They overlap
+** 2 P1 is completely contained within P2
+** 3 P2 is completely contained within P1
+** 4 P1 and P2 are the same polygon
+** NULL Either P1 or P2 or both are not valid polygons
+*/
+static void geopolyOverlapFunc(
+ sqlite3_context *context,
+ int argc,
+ sqlite3_value **argv
+){
+ GeoPoly *p1 = geopolyFuncParam(context, argv[0], 0);
+ GeoPoly *p2 = geopolyFuncParam(context, argv[1], 0);
+ if( p1 && p2 ){
+ int x = geopolyOverlap(p1, p2);
+ if( x<0 ){
+ sqlite3_result_error_nomem(context);
+ }else{
+ sqlite3_result_int(context, x);
+ }
+ }
+ sqlite3_free(p1);
+ sqlite3_free(p2);
+}
+
+/*
+** Enable or disable debugging output
+*/
+static void geopolyDebugFunc(
+ sqlite3_context *context,
+ int argc,
+ sqlite3_value **argv
+){
+#ifdef GEOPOLY_ENABLE_DEBUG
+ geo_debug = sqlite3_value_int(argv[0]);
+#endif
+}
+
+/*
+** This function is the implementation of both the xConnect and xCreate
+** methods of the geopoly virtual table.
+**
+** argv[0] -> module name
+** argv[1] -> database name
+** argv[2] -> table name
+** argv[...] -> column names...
+*/
+static int geopolyInit(
+ sqlite3 *db, /* Database connection */
+ void *pAux, /* One of the RTREE_COORD_* constants */
+ int argc, const char *const*argv, /* Parameters to CREATE TABLE statement */
+ sqlite3_vtab **ppVtab, /* OUT: New virtual table */
+ char **pzErr, /* OUT: Error message, if any */
+ int isCreate /* True for xCreate, false for xConnect */
+){
+ int rc = SQLITE_OK;
+ Rtree *pRtree;
+ sqlite3_int64 nDb; /* Length of string argv[1] */
+ sqlite3_int64 nName; /* Length of string argv[2] */
+ sqlite3_str *pSql;
+ char *zSql;
+ int ii;
+
+ sqlite3_vtab_config(db, SQLITE_VTAB_CONSTRAINT_SUPPORT, 1);
+
+ /* Allocate the sqlite3_vtab structure */
+ nDb = strlen(argv[1]);
+ nName = strlen(argv[2]);
+ pRtree = (Rtree *)sqlite3_malloc64(sizeof(Rtree)+nDb+nName+2);
+ if( !pRtree ){
+ return SQLITE_NOMEM;
+ }
+ memset(pRtree, 0, sizeof(Rtree)+nDb+nName+2);
+ pRtree->nBusy = 1;
+ pRtree->base.pModule = &rtreeModule;
+ pRtree->zDb = (char *)&pRtree[1];
+ pRtree->zName = &pRtree->zDb[nDb+1];
+ pRtree->eCoordType = RTREE_COORD_REAL32;
+ pRtree->nDim = 2;
+ pRtree->nDim2 = 4;
+ memcpy(pRtree->zDb, argv[1], nDb);
+ memcpy(pRtree->zName, argv[2], nName);
+
+
+ /* Create/Connect to the underlying relational database schema. If
+ ** that is successful, call sqlite3_declare_vtab() to configure
+ ** the r-tree table schema.
+ */
+ pSql = sqlite3_str_new(db);
+ sqlite3_str_appendf(pSql, "CREATE TABLE x(_shape");
+ pRtree->nAux = 1; /* Add one for _shape */
+ pRtree->nAuxNotNull = 1; /* The _shape column is always not-null */
+ for(ii=3; ii<argc; ii++){
+ pRtree->nAux++;
+ sqlite3_str_appendf(pSql, ",%s", argv[ii]);
+ }
+ sqlite3_str_appendf(pSql, ");");
+ zSql = sqlite3_str_finish(pSql);
+ if( !zSql ){
+ rc = SQLITE_NOMEM;
+ }else if( SQLITE_OK!=(rc = sqlite3_declare_vtab(db, zSql)) ){
+ *pzErr = sqlite3_mprintf("%s", sqlite3_errmsg(db));
+ }
+ sqlite3_free(zSql);
+ if( rc ) goto geopolyInit_fail;
+ pRtree->nBytesPerCell = 8 + pRtree->nDim2*4;
+
+ /* Figure out the node size to use. */
+ rc = getNodeSize(db, pRtree, isCreate, pzErr);
+ if( rc ) goto geopolyInit_fail;
+ rc = rtreeSqlInit(pRtree, db, argv[1], argv[2], isCreate);
+ if( rc ){
+ *pzErr = sqlite3_mprintf("%s", sqlite3_errmsg(db));
+ goto geopolyInit_fail;
+ }
+
+ *ppVtab = (sqlite3_vtab *)pRtree;
+ return SQLITE_OK;
+
+geopolyInit_fail:
+ if( rc==SQLITE_OK ) rc = SQLITE_ERROR;
+ assert( *ppVtab==0 );
+ assert( pRtree->nBusy==1 );
+ rtreeRelease(pRtree);
+ return rc;
+}
+
+
+/*
+** GEOPOLY virtual table module xCreate method.
+*/
+static int geopolyCreate(
+ sqlite3 *db,
+ void *pAux,
+ int argc, const char *const*argv,
+ sqlite3_vtab **ppVtab,
+ char **pzErr
+){
+ return geopolyInit(db, pAux, argc, argv, ppVtab, pzErr, 1);
+}
+
+/*
+** GEOPOLY virtual table module xConnect method.
+*/
+static int geopolyConnect(
+ sqlite3 *db,
+ void *pAux,
+ int argc, const char *const*argv,
+ sqlite3_vtab **ppVtab,
+ char **pzErr
+){
+ return geopolyInit(db, pAux, argc, argv, ppVtab, pzErr, 0);
+}
+
+
+/*
+** GEOPOLY virtual table module xFilter method.
+**
+** Query plans:
+**
+** 1 rowid lookup
+** 2 search for objects overlapping the same bounding box
+** that contains polygon argv[0]
+** 3 search for objects overlapping the same bounding box
+** that contains polygon argv[0]
+** 4 full table scan
+*/
+static int geopolyFilter(
+ sqlite3_vtab_cursor *pVtabCursor, /* The cursor to initialize */
+ int idxNum, /* Query plan */
+ const char *idxStr, /* Not Used */
+ int argc, sqlite3_value **argv /* Parameters to the query plan */
+){
+ Rtree *pRtree = (Rtree *)pVtabCursor->pVtab;
+ RtreeCursor *pCsr = (RtreeCursor *)pVtabCursor;
+ RtreeNode *pRoot = 0;
+ int rc = SQLITE_OK;
+ int iCell = 0;
+
+ rtreeReference(pRtree);
+
+ /* Reset the cursor to the same state as rtreeOpen() leaves it in. */
+ resetCursor(pCsr);
+
+ pCsr->iStrategy = idxNum;
+ if( idxNum==1 ){
+ /* Special case - lookup by rowid. */
+ RtreeNode *pLeaf; /* Leaf on which the required cell resides */
+ RtreeSearchPoint *p; /* Search point for the leaf */
+ i64 iRowid = sqlite3_value_int64(argv[0]);
+ i64 iNode = 0;
+ rc = findLeafNode(pRtree, iRowid, &pLeaf, &iNode);
+ if( rc==SQLITE_OK && pLeaf!=0 ){
+ p = rtreeSearchPointNew(pCsr, RTREE_ZERO, 0);
+ assert( p!=0 ); /* Always returns pCsr->sPoint */
+ pCsr->aNode[0] = pLeaf;
+ p->id = iNode;
+ p->eWithin = PARTLY_WITHIN;
+ rc = nodeRowidIndex(pRtree, pLeaf, iRowid, &iCell);
+ p->iCell = (u8)iCell;
+ RTREE_QUEUE_TRACE(pCsr, "PUSH-F1:");
+ }else{
+ pCsr->atEOF = 1;
+ }
+ }else{
+ /* Normal case - r-tree scan. Set up the RtreeCursor.aConstraint array
+ ** with the configured constraints.
+ */
+ rc = nodeAcquire(pRtree, 1, 0, &pRoot);
+ if( rc==SQLITE_OK && idxNum<=3 ){
+ RtreeCoord bbox[4];
+ RtreeConstraint *p;
+ assert( argc==1 );
+ assert( argv[0]!=0 );
+ geopolyBBox(0, argv[0], bbox, &rc);
+ if( rc ){
+ goto geopoly_filter_end;
+ }
+ pCsr->aConstraint = p = sqlite3_malloc(sizeof(RtreeConstraint)*4);
+ pCsr->nConstraint = 4;
+ if( p==0 ){
+ rc = SQLITE_NOMEM;
+ }else{
+ memset(pCsr->aConstraint, 0, sizeof(RtreeConstraint)*4);
+ memset(pCsr->anQueue, 0, sizeof(u32)*(pRtree->iDepth + 1));
+ if( idxNum==2 ){
+ /* Overlap query */
+ p->op = 'B';
+ p->iCoord = 0;
+ p->u.rValue = bbox[1].f;
+ p++;
+ p->op = 'D';
+ p->iCoord = 1;
+ p->u.rValue = bbox[0].f;
+ p++;
+ p->op = 'B';
+ p->iCoord = 2;
+ p->u.rValue = bbox[3].f;
+ p++;
+ p->op = 'D';
+ p->iCoord = 3;
+ p->u.rValue = bbox[2].f;
+ }else{
+ /* Within query */
+ p->op = 'D';
+ p->iCoord = 0;
+ p->u.rValue = bbox[0].f;
+ p++;
+ p->op = 'B';
+ p->iCoord = 1;
+ p->u.rValue = bbox[1].f;
+ p++;
+ p->op = 'D';
+ p->iCoord = 2;
+ p->u.rValue = bbox[2].f;
+ p++;
+ p->op = 'B';
+ p->iCoord = 3;
+ p->u.rValue = bbox[3].f;
+ }
+ }
+ }
+ if( rc==SQLITE_OK ){
+ RtreeSearchPoint *pNew;
+ pNew = rtreeSearchPointNew(pCsr, RTREE_ZERO, (u8)(pRtree->iDepth+1));
+ if( pNew==0 ){
+ rc = SQLITE_NOMEM;
+ goto geopoly_filter_end;
+ }
+ pNew->id = 1;
+ pNew->iCell = 0;
+ pNew->eWithin = PARTLY_WITHIN;
+ assert( pCsr->bPoint==1 );
+ pCsr->aNode[0] = pRoot;
+ pRoot = 0;
+ RTREE_QUEUE_TRACE(pCsr, "PUSH-Fm:");
+ rc = rtreeStepToLeaf(pCsr);
+ }
+ }
+
+geopoly_filter_end:
+ nodeRelease(pRtree, pRoot);
+ rtreeRelease(pRtree);
+ return rc;
+}
+
+/*
+** Rtree virtual table module xBestIndex method. There are three
+** table scan strategies to choose from (in order from most to
+** least desirable):
+**
+** idxNum idxStr Strategy
+** ------------------------------------------------
+** 1 "rowid" Direct lookup by rowid.
+** 2 "rtree" R-tree overlap query using geopoly_overlap()
+** 3 "rtree" R-tree within query using geopoly_within()
+** 4 "fullscan" full-table scan.
+** ------------------------------------------------
+*/
+static int geopolyBestIndex(sqlite3_vtab *tab, sqlite3_index_info *pIdxInfo){
+ int ii;
+ int iRowidTerm = -1;
+ int iFuncTerm = -1;
+ int idxNum = 0;
+
+ for(ii=0; ii<pIdxInfo->nConstraint; ii++){
+ struct sqlite3_index_constraint *p = &pIdxInfo->aConstraint[ii];
+ if( !p->usable ) continue;
+ if( p->iColumn<0 && p->op==SQLITE_INDEX_CONSTRAINT_EQ ){
+ iRowidTerm = ii;
+ break;
+ }
+ if( p->iColumn==0 && p->op>=SQLITE_INDEX_CONSTRAINT_FUNCTION ){
+ /* p->op==SQLITE_INDEX_CONSTRAINT_FUNCTION for geopoly_overlap()
+ ** p->op==(SQLITE_INDEX_CONTRAINT_FUNCTION+1) for geopoly_within().
+ ** See geopolyFindFunction() */
+ iFuncTerm = ii;
+ idxNum = p->op - SQLITE_INDEX_CONSTRAINT_FUNCTION + 2;
+ }
+ }
+
+ if( iRowidTerm>=0 ){
+ pIdxInfo->idxNum = 1;
+ pIdxInfo->idxStr = "rowid";
+ pIdxInfo->aConstraintUsage[iRowidTerm].argvIndex = 1;
+ pIdxInfo->aConstraintUsage[iRowidTerm].omit = 1;
+ pIdxInfo->estimatedCost = 30.0;
+ pIdxInfo->estimatedRows = 1;
+ pIdxInfo->idxFlags = SQLITE_INDEX_SCAN_UNIQUE;
+ return SQLITE_OK;
+ }
+ if( iFuncTerm>=0 ){
+ pIdxInfo->idxNum = idxNum;
+ pIdxInfo->idxStr = "rtree";
+ pIdxInfo->aConstraintUsage[iFuncTerm].argvIndex = 1;
+ pIdxInfo->aConstraintUsage[iFuncTerm].omit = 0;
+ pIdxInfo->estimatedCost = 300.0;
+ pIdxInfo->estimatedRows = 10;
+ return SQLITE_OK;
+ }
+ pIdxInfo->idxNum = 4;
+ pIdxInfo->idxStr = "fullscan";
+ pIdxInfo->estimatedCost = 3000000.0;
+ pIdxInfo->estimatedRows = 100000;
+ return SQLITE_OK;
+}
+
+
+/*
+** GEOPOLY virtual table module xColumn method.
+*/
+static int geopolyColumn(sqlite3_vtab_cursor *cur, sqlite3_context *ctx, int i){
+ Rtree *pRtree = (Rtree *)cur->pVtab;
+ RtreeCursor *pCsr = (RtreeCursor *)cur;
+ RtreeSearchPoint *p = rtreeSearchPointFirst(pCsr);
+ int rc = SQLITE_OK;
+ RtreeNode *pNode = rtreeNodeOfFirstSearchPoint(pCsr, &rc);
+
+ if( rc ) return rc;
+ if( p==0 ) return SQLITE_OK;
+ if( i==0 && sqlite3_vtab_nochange(ctx) ) return SQLITE_OK;
+ if( i<=pRtree->nAux ){
+ if( !pCsr->bAuxValid ){
+ if( pCsr->pReadAux==0 ){
+ rc = sqlite3_prepare_v3(pRtree->db, pRtree->zReadAuxSql, -1, 0,
+ &pCsr->pReadAux, 0);
+ if( rc ) return rc;
+ }
+ sqlite3_bind_int64(pCsr->pReadAux, 1,
+ nodeGetRowid(pRtree, pNode, p->iCell));
+ rc = sqlite3_step(pCsr->pReadAux);
+ if( rc==SQLITE_ROW ){
+ pCsr->bAuxValid = 1;
+ }else{
+ sqlite3_reset(pCsr->pReadAux);
+ if( rc==SQLITE_DONE ) rc = SQLITE_OK;
+ return rc;
+ }
+ }
+ sqlite3_result_value(ctx, sqlite3_column_value(pCsr->pReadAux, i+2));
+ }
+ return SQLITE_OK;
+}
+
+
+/*
+** The xUpdate method for GEOPOLY module virtual tables.
+**
+** For DELETE:
+**
+** argv[0] = the rowid to be deleted
+**
+** For INSERT:
+**
+** argv[0] = SQL NULL
+** argv[1] = rowid to insert, or an SQL NULL to select automatically
+** argv[2] = _shape column
+** argv[3] = first application-defined column....
+**
+** For UPDATE:
+**
+** argv[0] = rowid to modify. Never NULL
+** argv[1] = rowid after the change. Never NULL
+** argv[2] = new value for _shape
+** argv[3] = new value for first application-defined column....
+*/
+static int geopolyUpdate(
+ sqlite3_vtab *pVtab,
+ int nData,
+ sqlite3_value **aData,
+ sqlite_int64 *pRowid
+){
+ Rtree *pRtree = (Rtree *)pVtab;
+ int rc = SQLITE_OK;
+ RtreeCell cell; /* New cell to insert if nData>1 */
+ i64 oldRowid; /* The old rowid */
+ int oldRowidValid; /* True if oldRowid is valid */
+ i64 newRowid; /* The new rowid */
+ int newRowidValid; /* True if newRowid is valid */
+ int coordChange = 0; /* Change in coordinates */
+
+ if( pRtree->nNodeRef ){
+ /* Unable to write to the btree while another cursor is reading from it,
+ ** since the write might do a rebalance which would disrupt the read
+ ** cursor. */
+ return SQLITE_LOCKED_VTAB;
+ }
+ rtreeReference(pRtree);
+ assert(nData>=1);
+
+ oldRowidValid = sqlite3_value_type(aData[0])!=SQLITE_NULL;;
+ oldRowid = oldRowidValid ? sqlite3_value_int64(aData[0]) : 0;
+ newRowidValid = nData>1 && sqlite3_value_type(aData[1])!=SQLITE_NULL;
+ newRowid = newRowidValid ? sqlite3_value_int64(aData[1]) : 0;
+ cell.iRowid = newRowid;
+
+ if( nData>1 /* not a DELETE */
+ && (!oldRowidValid /* INSERT */
+ || !sqlite3_value_nochange(aData[2]) /* UPDATE _shape */
+ || oldRowid!=newRowid) /* Rowid change */
+ ){
+ assert( aData[2]!=0 );
+ geopolyBBox(0, aData[2], cell.aCoord, &rc);
+ if( rc ){
+ if( rc==SQLITE_ERROR ){
+ pVtab->zErrMsg =
+ sqlite3_mprintf("_shape does not contain a valid polygon");
+ }
+ goto geopoly_update_end;
+ }
+ coordChange = 1;
+
+ /* If a rowid value was supplied, check if it is already present in
+ ** the table. If so, the constraint has failed. */
+ if( newRowidValid && (!oldRowidValid || oldRowid!=newRowid) ){
+ int steprc;
+ sqlite3_bind_int64(pRtree->pReadRowid, 1, cell.iRowid);
+ steprc = sqlite3_step(pRtree->pReadRowid);
+ rc = sqlite3_reset(pRtree->pReadRowid);
+ if( SQLITE_ROW==steprc ){
+ if( sqlite3_vtab_on_conflict(pRtree->db)==SQLITE_REPLACE ){
+ rc = rtreeDeleteRowid(pRtree, cell.iRowid);
+ }else{
+ rc = rtreeConstraintError(pRtree, 0);
+ }
+ }
+ }
+ }
+
+ /* If aData[0] is not an SQL NULL value, it is the rowid of a
+ ** record to delete from the r-tree table. The following block does
+ ** just that.
+ */
+ if( rc==SQLITE_OK && (nData==1 || (coordChange && oldRowidValid)) ){
+ rc = rtreeDeleteRowid(pRtree, oldRowid);
+ }
+
+ /* If the aData[] array contains more than one element, elements
+ ** (aData[2]..aData[argc-1]) contain a new record to insert into
+ ** the r-tree structure.
+ */
+ if( rc==SQLITE_OK && nData>1 && coordChange ){
+ /* Insert the new record into the r-tree */
+ RtreeNode *pLeaf = 0;
+ if( !newRowidValid ){
+ rc = rtreeNewRowid(pRtree, &cell.iRowid);
+ }
+ *pRowid = cell.iRowid;
+ if( rc==SQLITE_OK ){
+ rc = ChooseLeaf(pRtree, &cell, 0, &pLeaf);
+ }
+ if( rc==SQLITE_OK ){
+ int rc2;
+ pRtree->iReinsertHeight = -1;
+ rc = rtreeInsertCell(pRtree, pLeaf, &cell, 0);
+ rc2 = nodeRelease(pRtree, pLeaf);
+ if( rc==SQLITE_OK ){
+ rc = rc2;
+ }
+ }
+ }
+
+ /* Change the data */
+ if( rc==SQLITE_OK && nData>1 ){
+ sqlite3_stmt *pUp = pRtree->pWriteAux;
+ int jj;
+ int nChange = 0;
+ sqlite3_bind_int64(pUp, 1, cell.iRowid);
+ assert( pRtree->nAux>=1 );
+ if( sqlite3_value_nochange(aData[2]) ){
+ sqlite3_bind_null(pUp, 2);
+ }else{
+ GeoPoly *p = 0;
+ if( sqlite3_value_type(aData[2])==SQLITE_TEXT
+ && (p = geopolyFuncParam(0, aData[2], &rc))!=0
+ && rc==SQLITE_OK
+ ){
+ sqlite3_bind_blob(pUp, 2, p->hdr, 4+8*p->nVertex, SQLITE_TRANSIENT);
+ }else{
+ sqlite3_bind_value(pUp, 2, aData[2]);
+ }
+ sqlite3_free(p);
+ nChange = 1;
+ }
+ for(jj=1; jj<nData-2; jj++){
+ nChange++;
+ sqlite3_bind_value(pUp, jj+2, aData[jj+2]);
+ }
+ if( nChange ){
+ sqlite3_step(pUp);
+ rc = sqlite3_reset(pUp);
+ }
+ }
+
+geopoly_update_end:
+ rtreeRelease(pRtree);
+ return rc;
+}
+
+/*
+** Report that geopoly_overlap() is an overloaded function suitable
+** for use in xBestIndex.
+*/
+static int geopolyFindFunction(
+ sqlite3_vtab *pVtab,
+ int nArg,
+ const char *zName,
+ void (**pxFunc)(sqlite3_context*,int,sqlite3_value**),
+ void **ppArg
+){
+ if( sqlite3_stricmp(zName, "geopoly_overlap")==0 ){
+ *pxFunc = geopolyOverlapFunc;
+ *ppArg = 0;
+ return SQLITE_INDEX_CONSTRAINT_FUNCTION;
+ }
+ if( sqlite3_stricmp(zName, "geopoly_within")==0 ){
+ *pxFunc = geopolyWithinFunc;
+ *ppArg = 0;
+ return SQLITE_INDEX_CONSTRAINT_FUNCTION+1;
+ }
+ return 0;
+}
+
+
+static sqlite3_module geopolyModule = {
+ 3, /* iVersion */
+ geopolyCreate, /* xCreate - create a table */
+ geopolyConnect, /* xConnect - connect to an existing table */
+ geopolyBestIndex, /* xBestIndex - Determine search strategy */
+ rtreeDisconnect, /* xDisconnect - Disconnect from a table */
+ rtreeDestroy, /* xDestroy - Drop a table */
+ rtreeOpen, /* xOpen - open a cursor */
+ rtreeClose, /* xClose - close a cursor */
+ geopolyFilter, /* xFilter - configure scan constraints */
+ rtreeNext, /* xNext - advance a cursor */
+ rtreeEof, /* xEof */
+ geopolyColumn, /* xColumn - read data */
+ rtreeRowid, /* xRowid - read data */
+ geopolyUpdate, /* xUpdate - write data */
+ rtreeBeginTransaction, /* xBegin - begin transaction */
+ rtreeEndTransaction, /* xSync - sync transaction */
+ rtreeEndTransaction, /* xCommit - commit transaction */
+ rtreeEndTransaction, /* xRollback - rollback transaction */
+ geopolyFindFunction, /* xFindFunction - function overloading */
+ rtreeRename, /* xRename - rename the table */
+ rtreeSavepoint, /* xSavepoint */
+ 0, /* xRelease */
+ 0, /* xRollbackTo */
+ rtreeShadowName /* xShadowName */
+};
+
+static int sqlite3_geopoly_init(sqlite3 *db){
+ int rc = SQLITE_OK;
+ static const struct {
+ void (*xFunc)(sqlite3_context*,int,sqlite3_value**);
+ signed char nArg;
+ unsigned char bPure;
+ const char *zName;
+ } aFunc[] = {
+ { geopolyAreaFunc, 1, 1, "geopoly_area" },
+ { geopolyBlobFunc, 1, 1, "geopoly_blob" },
+ { geopolyJsonFunc, 1, 1, "geopoly_json" },
+ { geopolySvgFunc, -1, 1, "geopoly_svg" },
+ { geopolyWithinFunc, 2, 1, "geopoly_within" },
+ { geopolyContainsPointFunc, 3, 1, "geopoly_contains_point" },
+ { geopolyOverlapFunc, 2, 1, "geopoly_overlap" },
+ { geopolyDebugFunc, 1, 0, "geopoly_debug" },
+ { geopolyBBoxFunc, 1, 1, "geopoly_bbox" },
+ { geopolyXformFunc, 7, 1, "geopoly_xform" },
+ { geopolyRegularFunc, 4, 1, "geopoly_regular" },
+ { geopolyCcwFunc, 1, 1, "geopoly_ccw" },
+ };
+ static const struct {
+ void (*xStep)(sqlite3_context*,int,sqlite3_value**);
+ void (*xFinal)(sqlite3_context*);
+ const char *zName;
+ } aAgg[] = {
+ { geopolyBBoxStep, geopolyBBoxFinal, "geopoly_group_bbox" },
+ };
+ int i;
+ for(i=0; i<sizeof(aFunc)/sizeof(aFunc[0]) && rc==SQLITE_OK; i++){
+ int enc;
+ if( aFunc[i].bPure ){
+ enc = SQLITE_UTF8|SQLITE_DETERMINISTIC|SQLITE_INNOCUOUS;
+ }else{
+ enc = SQLITE_UTF8|SQLITE_DIRECTONLY;
+ }
+ rc = sqlite3_create_function(db, aFunc[i].zName, aFunc[i].nArg,
+ enc, 0,
+ aFunc[i].xFunc, 0, 0);
+ }
+ for(i=0; i<sizeof(aAgg)/sizeof(aAgg[0]) && rc==SQLITE_OK; i++){
+ rc = sqlite3_create_function(db, aAgg[i].zName, 1,
+ SQLITE_UTF8|SQLITE_DETERMINISTIC|SQLITE_INNOCUOUS, 0,
+ 0, aAgg[i].xStep, aAgg[i].xFinal);
+ }
+ if( rc==SQLITE_OK ){
+ rc = sqlite3_create_module_v2(db, "geopoly", &geopolyModule, 0, 0);
+ }
+ return rc;
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-12-24 13:58:44 +0000