mirrors/sqlite
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
e9f740e484
sqlite/ext/rtree/geopoly.c
drh f98109e7cb In the xUpdate method of the GeoPoly virtual table, make sure that the 
number of updated columns does not exceed the underlying implementation, even
if the virtual table object records an excess number of column in the nAux
field due to table constraints in the table definition.
Fix for the problem reported by
[forum:/forumpost/a096ab7d96bb057a|forum post a096ab7d96bb057a].

FossilOrigin-Name: a04dd475c2c8afdded0afecfd34c6c4c2f48cce520e09b7fa5172ff30f09209e
2022-08-01 13:14:03 +00:00

1815 lines
51 KiB
C
Raw Blame History

/*
** 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;
}
Reference in New Issue View Git Blame Copy Permalink