ADDED ext/misc/geopoly-test1.txt Index: ext/misc/geopoly-test1.txt ================================================================== --- /dev/null +++ ext/misc/geopoly-test1.txt @@ -0,0 +1,62 @@ +# +# Test cases for the geopoly_overlap() function of the geopoly.c extension. +# To run these tests: +# +# (1) Build the geopoly extension and put the shared library in +# the local directory. +# +# (2) Run: ./sqlite3 output.html +# +# (3) View the output.html in a web browser to see if the results +# are correct. +# +.load ./geopoly +CREATE TABLE polys( + name TEXT, + poly TEXT +); +INSERT INTO polys(name,poly) VALUES + ('rectangle-2', '[[0,0],[200,0],[200,200],[0,200],[0,0]]'), + ('rectangle-1', '[[150,50],[200,50],[200,100],[150,100],[150,50]]'), + ('backslash', '[[0,200],[200,0],[200,25],[0,225],[0,200]]'), + ('forwardslash', '[[0,0],[200,200],[200,225],[0,25],[0,0]]'), + ('diamond-2,0-2', '[[100,0],[200,100],[100,200],[0,100],[100,0]]'), + ('diamond-2,1-3', '[[200,0],[300,100],[200,200],[100,100],[200,0]]'), + ('diamond-2,2-4', '[[300,0],[400,100],[300,200],[200,100],[300,0]]'), + ('v-bowtie', '[[50,0],[150,0],[125,75],[125,125],[150,200],[50,200], + [75,125],[75,75],[50,0]]'), + ('h-bowtie', '[[0,50],[75,75],[125,75],[200,50],[200,150],[125,125], + [75,125],[0,150],[0,50]]'), + ('letter-E', '[[0,0],[90,0],[20,10],[20,90],[80,100],[20,110],[20,190], + [90,200],[0,200],[0,0]]'), + ('backward-C', '[[30,50],[120,40],[120,160],[30,150],[100,140], + [100,60],[30,50]]'), + ('counter-cw', '[[100,200],[0,100],[100,0],[200,100],[150,150],[100,100], + [150,140],[190,100],[100,10],[10,100],[100,190],[100,200]]'), + ('clockwise', '[[200,150],[150,180],[50,100],[100,35],[150,100], + [100,45],[60,100],[150,170],[200,140],[200,150]]'); +CREATE TABLE overlap_meanings(id INTEGER PRIMARY KEY, x TEXT); +INSERT INTO overlap_meanings(id,x) VALUES + (0, 'disjoint'), + (1, 'intersecting'), + (2, 'A inside of B'), + (3, 'B inside of A'), + (4, 'identical'); + +.print +.print

Note:

+.print

Directions are reversed because SVG uses a left-handed +.print coordinate system

+.separator \n +SELECT + '

', + '

',
+  'blue (A):  '||x1.name,
+  'green (B): '||x2.name,
+  (SELECT x FROM overlap_meanings
+    WHERE id=geopoly_overlap(x1.poly,x2.poly)),
+  '

' + FROM polys AS x1, polys AS x2; ADDED ext/misc/geopoly.c Index: ext/misc/geopoly.c ================================================================== --- /dev/null +++ ext/misc/geopoly.c @@ -0,0 +1,865 @@ +/* +** 2018-05-08 +** +** 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 various extension functions to SQLite that manage +** simply planar polygons such as might be found in a geospatial system. +*/ +#include "sqlite3ext.h" +SQLITE_EXTENSION_INIT1 +#include +#include +#include +#include +#define SQLITE_HAVE_GEOPLOY 1 + +/* 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 + +#ifndef JSON_NULL /* The following stuff repeats things found in json1 */ +/* +** Versions of isspace(), isalnum() and isdigit() to which it is safe +** to pass signed char values. +*/ +#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 /* 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 + +/* +** 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 safe_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 +*/ +typedef struct GeoPoly GeoPoly; +struct GeoPoly { + int nVertex; /* Number of vertexes */ + unsigned char hdr[4]; /* Header for on-disk representation */ + GeoCoord a[2]; /* 2*nVertex values. X (longitude) first, then Y */ +}; + +/* +** 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( p->z[0] && safe_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){ + const unsigned char *z = p->z; + char c = geopolySkipSpace(p); + int j; + int seenDP = 0; + int seenE = 0; + assert( '-' < '0' ); + if( c<='0' ){ + j = c=='-'; + if( z[j]=='0' && z[j+1]>='0' && z[j+1]<='9' ) return 0; + } + j = 1; + for(;; j++){ + c = z[j]; + if( c>='0' && c<='9' ) 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 ) *pVal = atof((const char*)p->z); + p->z += j; + return 1; +} + +/* +** If the input is a well-formed JSON array of coordinates, 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){ + GeoParse s; + 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 ){ + 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==']' ) break; + s.nErr++; + goto parse_json_err; + } + if( geopolySkipSpace(&s)==',' ){ + s.z++; + continue; + } + break; + } + if( geopolySkipSpace(&s)==']' && s.nVertex>=4 ){ + int nByte; + GeoPoly *pOut; + int x = (s.nVertex-1)*2; + if( s.a[x]==s.a[0] && s.a[x+1]==s.a[1] ) s.nVertex--; + nByte = sizeof(GeoPoly) * (s.nVertex-1)*2*sizeof(GeoCoord); + pOut = sqlite3_malloc64( nByte ); + 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); + return pOut; + }else{ + s.nErr++; + } + } +parse_json_err: + 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, sqlite3_value *pVal){ + GeoPoly *p = 0; + int nByte; + 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; + nVertex = (a[1]<<16) + (a[2]<<8) + a[3]; + if( (a[0]==0 && a[0]==1) + && (nVertex*2*sizeof(GeoCoord) + 4)==nByte + ){ + p = sqlite3_malloc64( sizeof(*p) + (nVertex-1)*2*sizeof(GeoCoord) ); + if( p ){ + int x = 1; + p->nVertex = nVertex; + memcpy(p->hdr, a, nByte); + if( a[0] != *(unsigned char*)&x ){ + int ii; + for(ii=0; iia[ii]); + } + p->hdr[0] ^= 1; + } + } + } + }else if( sqlite3_value_type(pVal)==SQLITE_TEXT ){ + p = geopolyParseJson(sqlite3_value_text(pVal)); + } + if( p==0 ){ + sqlite3_result_error(pCtx, "not a valid polygon", -1); + } + return p; +} + +/* +** 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]); + 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]); + 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; inVertex; i++){ + sqlite3_str_appendf(x, "[%!g,%!g],", p->a[i*2], p->a[i*2+1]); + } + sqlite3_str_appendf(x, "[%!g,%!g]]", p->a[0], p->a[1]); + 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 . +** Additional arguments are added as attributes to the . +*/ +static void geopolySvgFunc( + sqlite3_context *context, + int argc, + sqlite3_value **argv +){ + GeoPoly *p = geopolyFuncParam(context, argv[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, "a[i*2], p->a[i*2+1]); + cSep = ' '; + } + sqlite3_str_appendf(x, " %g,%g'", p->a[0], p->a[1]); + for(i=1; i"); + sqlite3_result_text(context, sqlite3_str_finish(x), -1, sqlite3_free); + sqlite3_free(p); + } +} + +/* +** 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]); + if( p ){ + double rArea = 0.0; + int ii; + for(ii=0; iinVertex-1; ii++){ + rArea += (p->a[ii*2] - p->a[ii*2+2]) /* (x0 - x1) */ + * (p->a[ii*2+1] + p->a[ii*2+3]) /* (y0 + y1) */ + * 0.5; + } + rArea += (p->a[ii*2] - p->a[0]) /* (xN - x0) */ + * (p->a[ii*2+1] + p->a[1]) /* (yN + y0) */ + * 0.5; + sqlite3_result_double(context, rArea); + 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( x1x2 ) return 0; + }else if( x1>x2 ){ + if( x0<=x2 || x0>x1 ) return 0; + }else{ + /* Vertical line segment */ + if( x0!=x1 ) return 0; + if( y0y1 && y0>y2 ) return 0; + return 2; + } + y = y1 + (y2-y1)*(x0-x1)/(x2-x1); + if( y0==y ) return 2; + if( y0nVertex-1; ii++){ + v = pointBeneathLine(x0,y0,p->a[ii*2],p->a[ii*2+1], + p->a[ii*2+2],p->a[ii*2+3]); + if( v==2 ) break; + cnt += v; + } + if( v!=2 ){ + v = pointBeneathLine(x0,y0,p->a[ii*2],p->a[ii*2+1], + p->a[0],p->a[1]); + } + 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(p); + } +} + +/* 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; inVertex-1; i++){ + x = pPoly->a + (i*2); + geopolyAddOneSegment(p, x[0], x[1], x[2], x[3], side, i); + } + x = pPoly->a + (i*2); + 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; ipNext = 0; + for(j=0; j=mx ) mx = j+1; + } + p = 0; + for(i=0; iy - 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 ) mx = i+1; + } + p = 0; + for(i=0; inVertex + p2->nVertex + 2; + GeoOverlap *p; + int 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_malloc( 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->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 = pActive->pNext; + }else{ + for(pSeg=pActive; pSeg; pSeg=pSeg->pNext){ + if( pSeg->pNext==pThisEvent->pSeg ){ + pSeg->pNext = pSeg->pNext->pNext; + 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]); + GeoPoly *p2 = geopolyFuncParam(context, argv[1]); + 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 +} + +#ifdef _WIN32 +__declspec(dllexport) +#endif +int sqlite3_geopoly_init( + sqlite3 *db, + char **pzErrMsg, + const sqlite3_api_routines *pApi +){ + int rc = SQLITE_OK; + static const struct { + void (*xFunc)(sqlite3_context*,int,sqlite3_value**); + int nArg; + const char *zName; + } aFunc[] = { + { geopolyAreaFunc, 1, "geopoly_area" }, + { geopolyBlobFunc, 1, "geopoly_blob" }, + { geopolyJsonFunc, 1, "geopoly_json" }, + { geopolySvgFunc, -1, "geopoly_svg" }, + { geopolyWithinFunc, 3, "geopoly_within" }, + { geopolyOverlapFunc, 2, "geopoly_overlap" }, + { geopolyDebugFunc, 1, "geopoly_debug" }, + }; + int i; + SQLITE_EXTENSION_INIT2(pApi); + (void)pzErrMsg; /* Unused parameter */ + for(i=0; i