| /* |
| ** 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))>=(int)(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); |
| (void)argc; |
| 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); |
| (void)argc; |
| 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; |
| (void)argc; |
| 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); |
| (void)argc; |
| 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); |
| (void)argc; |
| 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; |
| (void)argc; |
| |
| 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*)ⅈ |
| 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); |
| (void)argc; |
| 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)argc; |
| (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; |
| (void)argc; |
| |
| 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); |
| (void)argc; |
| 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); |
| (void)argc; |
| 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 |
| ){ |
| (void)context; |
| (void)argc; |
| #ifdef GEOPOLY_ENABLE_DEBUG |
| geo_debug = sqlite3_value_int(argv[0]); |
| #else |
| (void)argv; |
| #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; |
| (void)pAux; |
| |
| sqlite3_vtab_config(db, SQLITE_VTAB_CONSTRAINT_SUPPORT, 1); |
| sqlite3_vtab_config(db, SQLITE_VTAB_INNOCUOUS); |
| |
| /* Allocate the sqlite3_vtab structure */ |
| nDb = strlen(argv[1]); |
| nName = strlen(argv[2]); |
| pRtree = (Rtree *)sqlite3_malloc64(sizeof(Rtree)+nDb+nName*2+8); |
| if( !pRtree ){ |
| return SQLITE_NOMEM; |
| } |
| memset(pRtree, 0, sizeof(Rtree)+nDb+nName*2+8); |
| pRtree->nBusy = 1; |
| pRtree->base.pModule = &rtreeModule; |
| pRtree->zDb = (char *)&pRtree[1]; |
| pRtree->zName = &pRtree->zDb[nDb+1]; |
| pRtree->zNodeName = &pRtree->zName[nName+1]; |
| pRtree->eCoordType = RTREE_COORD_REAL32; |
| pRtree->nDim = 2; |
| pRtree->nDim2 = 4; |
| memcpy(pRtree->zDb, argv[1], nDb); |
| memcpy(pRtree->zName, argv[2], nName); |
| memcpy(pRtree->zNodeName, argv[2], nName); |
| memcpy(&pRtree->zNodeName[nName], "_node", 6); |
| |
| |
| /* 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; |
| (void)idxStr; |
| |
| 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; |
| (void)tab; |
| |
| 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; |
| 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 |
| ){ |
| (void)pVtab; |
| (void)nArg; |
| 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 */ |
| rtreeIntegrity /* xIntegrity */ |
| }; |
| |
| 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" }, |
| }; |
| unsigned 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; |
| } |