Changeset: 84c2d53795b9 for MonetDB
URL: http://dev.monetdb.org/hg/MonetDB?cmd=changeset;node=84c2d53795b9
Modified Files:
geom/monetdb5/geom.h
geom/monetdb5/geom.mal
geom/monetdb5/geomPoints.c
geom/sql/40_geom.sql
monetdb5/optimizer/opt_geospatial.c
Branch: geo
Log Message:
implementation of distance for x,y coordinates without the use of geos
diffs (259 lines):
diff --git a/geom/monetdb5/geom.h b/geom/monetdb5/geom.h
--- a/geom/monetdb5/geom.h
+++ b/geom/monetdb5/geom.h
@@ -266,7 +266,8 @@ geom_export str wkbCoordinateFromMBR_bat
geom_export str wkbPointsContains1_geom_bat(bat* outBAT_id, wkb** geomWKB,
bat* xBAT_id, bat* yBAT_id, int* srid);
geom_export str wkbPointsContains2_geom_bat(bat* outBAT_id, wkb** geomWKB,
bat* xBAT_id, bat* yBAT_id, int* srid);
-geom_export str wkbPointsDistance_geom_bat(bat* outBAT_id, wkb** geomWKB, bat*
xBAT_id, bat* yBAT_id, int* srid);
+geom_export str wkbPointsDistance1_geom_bat(bat* outBAT_id, wkb** geomWKB,
bat* xBAT_id, bat* yBAT_id, int* srid);
+geom_export str wkbPointsDistance2_geom_bat(bat* outBAT_id, wkb** geomWKB,
bat* xBAT_id, bat* yBAT_id, int* srid);
geom_export str wkbFilterWithImprints_geom_bat(bat*, wkb**, bat*, bat*);
diff --git a/geom/monetdb5/geom.mal b/geom/monetdb5/geom.mal
--- a/geom/monetdb5/geom.mal
+++ b/geom/monetdb5/geom.mal
@@ -443,12 +443,13 @@ geom.prelude();
module batgeom;
+#the 1 version of the functions use geos the 2 versions have custom
implementations
+
command Distance(a:bat[:oid,:wkb], b:bat[:oid,:wkb]) :bat[:oid,:dbl] address
wkbDistance_bat;
command Distance(a:wkb, b:bat[:oid,:wkb]) :bat[:oid,:dbl] address
wkbDistance_geom_bat;
command Distance(a:bat[:oid,:wkb], b:wkb) :bat[:oid,:dbl] address
wkbDistance_bat_geom;
-command Distance(g:wkb, x:bat[:oid,:dbl], y:bat[:oid,:dbl], srid:int)
:bat[:oid,:dbl] address wkbPointsDistance_geom_bat;
-command Distance(g:wkb, x:bat[:oid,:dbl], y:bat[:oid,:dbl], o:bat[:oid,:oid],
srid:int) :bat[:oid,:dbl] address wkbFilteredPointsDistance_geom_bat;
-
+command Distance1(g:wkb, x:bat[:oid,:dbl], y:bat[:oid,:dbl], srid:int)
:bat[:oid,:dbl] address wkbPointsDistance1_geom_bat;
+command Distance2(g:wkb, x:bat[:oid,:dbl], y:bat[:oid,:dbl], srid:int)
:bat[:oid,:dbl] address wkbPointsDistance2_geom_bat;
command Contains(a:bat[:oid,:wkb], b:bat[:oid,:wkb]) :bat[:oid,:bit] address
wkbContains_bat;
command Contains(a:wkb, b:bat[:oid,:wkb]) :bat[:oid,:bit] address
wkbContains_geom_bat;
diff --git a/geom/monetdb5/geomPoints.c b/geom/monetdb5/geomPoints.c
--- a/geom/monetdb5/geomPoints.c
+++ b/geom/monetdb5/geomPoints.c
@@ -255,7 +255,6 @@ inline int isLeft( double P0x, double P0
#define isRight(x0, y0, x1, y1, x2, y2) isLeft(x0, y0, x1, y1, x2, y2)-1
-//static str pnpoly_(int *out, int nvert, dbl *vx, dbl *vy, int *point_x, int
*point_y) {
static str pnpoly_(int *out, const GEOSGeometry *geosGeometry, int *point_x,
int *point_y) {
BAT *bo = NULL, *bpx = NULL, *bpy = NULL;
dbl *px = NULL, *py = NULL;
@@ -432,16 +431,6 @@ static str pnpolyWithHoles_(int *out, GE
else if (yCurrent >= py[i] && yNext <= py[i]) {
wn+=isRight( xCurrent, yCurrent, xNext, yNext,
px[i], py[i]);
}
-
- // if (yCurrent <= py[i]) {
- // if (yNext > py[i])
- // if (isLeft( xCurrent, yCurrent, xNext,
yNext, px[i], py[i]) > 0)
- // ++wn;
- // } else {
- // if (yNext <= py[i])
- // if (isLeft(xCurrent, yCurrent, xNext,
yNext, px[i], py[i]) < 0)
- // --wn;
- // }
}
//It is in one of the holes no reason to check the others
@@ -468,15 +457,6 @@ static str pnpolyWithHoles_(int *out, GE
else if (yCurrent >= py[i] && yNext <= py[i]) {
wn+=isRight( xCurrent, yCurrent, xNext, yNext, px[i],
py[i]);
}
- //if (yCurrent <= py[i]) {
- // if (yNext > py[i])
- // if (isLeft( xCurrent, yCurrent, xNext,
yNext, px[i], py[i]) > 0)
- /// ++wn;
- //} else {
- // if (yNext <= py[i])
- // if (isLeft(xCurrent, yCurrent, xNext,
yNext, px[i], py[i]) < 0)
- // --wn;
- //}
}
*cs++ = wn&1;
}
@@ -541,7 +521,7 @@ static BAT* BATDistance(wkb** geomWKB, B
return NULL;
}
- //set the first idx of the new BAT equal to that of the x BAT (which is
equal to the y BAT)
+ //set the first idx of the new BAT equal to that of the geometries BAT
BATseqbase(outBAT, geometriesBAT->hseqbase);
//iterator over the BATs
@@ -564,7 +544,7 @@ static BAT* BATDistance(wkb** geomWKB, B
return outBAT;
}
-str wkbPointsDistance_geom_bat(bat* outBAT_id, wkb** geomWKB, bat* xBAT_id,
bat* yBAT_id, int* srid) {
+str wkbPointsDistance1_geom_bat(bat* outBAT_id, wkb** geomWKB, bat* xBAT_id,
bat* yBAT_id, int* srid) {
BAT *xBAT=NULL, *yBAT=NULL, *outBAT=NULL;
BAT *pointsBAT = NULL, *pointsWithSRIDBAT=NULL;
str ret=MAL_SUCCEED;
@@ -620,6 +600,123 @@ clean:
return ret;
}
+
+/* Alternative implementation of distance that computes the euclidean distance
of two points
+ * (other geometries are not yet supported)
+ * Using the Euclidean distance is appropriate when the geometries are
expressed in projected
+ * reference system, but what happens with the rest? Maybe we should check nd
report an error
+ * when a geographic srid is used. PostGIS, however, does not perform any
check on the srid
+ * and always computes the Euclidean distance
+ * */
+
+static BAT* point2point_distance(GEOSGeom geosGeometry, BAT *xBAT, BAT *yBAT) {
+ BAT *outBAT = NULL;
+ const GEOSCoordSequence *coordSeq;
+ double xPointCoordinate = 0.0, yPointCoordinate = 0.0;
+ double *xBATCoordinate = NULL, *yBATCoordinate = NULL, *distancesArray
= NULL;
+ BUN i=0;
+
+ //create a new BAT
+ if ((outBAT = BATnew(TYPE_void, ATOMindex("dbl"), BATcount(xBAT),
TRANSIENT)) == NULL) {
+ GDKerror("BATDistance: Could not create new BAT for the
output");
+ return NULL;
+ }
+
+ //set the first idx of the new BAT equal to that of the x BAT (which is
equal to the y BAT)
+ BATseqbase(outBAT, xBAT->hseqbase);
+
+
+ //the geometry is a point. Get the x, y coordinates of it
+ if(!(coordSeq = GEOSGeom_getCoordSeq(geosGeometry))) {
+ GDKerror("batgeom.Distance: GEOSGeom_getCoordSeq failed");
+ return NULL;
+ }
+ if(GEOSCoordSeq_getX(coordSeq, 0, &xPointCoordinate) == -1) {
+ GDKerror("batgeom.Distance: GEOSCoordSeq_getX failed");
+ return NULL;
+ }
+ if(GEOSCoordSeq_getY(coordSeq, 0, &yPointCoordinate) == -1) {
+ GDKerror("batgeom.Distance: GEOSCoordSeq_getY failed");
+ return NULL;
+ }
+
+ //get the x and y coordinates from the BATs (fixed size)
+ xBATCoordinate = (double*) Tloc(xBAT, BUNfirst(xBAT));
+ yBATCoordinate = (double*) Tloc(yBAT, BUNfirst(yBAT));
+ distancesArray = (double*) Tloc(outBAT,BUNfirst(outBAT));
+ //iterate
+ for (i = 0; i < BATcount(xBAT); i++) {
+ distancesArray[i] =
sqrt(pow((xBATCoordinate[i]-xPointCoordinate),
2.0)+pow((yBATCoordinate[i]-yPointCoordinate), 2.0));
+/* if(i%1000 == 0) {
+ fprintf(stderr, "%u: %f\n", (unsigned int)i,
distancesArray[i]);
+ fprintf(stderr, "\t(%f, %f), (%f, %f)\n",
xPointCoordinate, yPointCoordinate, xBATCoordinate[i], yBATCoordinate[i]);
+ }
+*/
+ }
+
+ BATsetcount(outBAT,BATcount(xBAT));
+ BATderiveProps(outBAT,FALSE);
+
+ return outBAT;
+
+
+}
+
+str wkbPointsDistance2_geom_bat(bat* outBAT_id, wkb** geomWKB, bat* xBAT_id,
bat* yBAT_id, int* srid) {
+ BAT *xBAT=NULL, *yBAT=NULL, *outBAT=NULL;
+ GEOSGeom geosGeometry;
+ str ret=MAL_SUCCEED;
+
+ //check if geometry a and the points have the same srid
+ if((*geomWKB)->srid != *srid)
+ return createException(MAL, "batgeom.Distance", "Geometry and
points should have the same srid");
+
+ //get the descriptors of the BATs
+ if ((xBAT = BATdescriptor(*xBAT_id)) == NULL) {
+ throw(MAL, "batgeom.Distance", RUNTIME_OBJECT_MISSING);
+ }
+ if ((yBAT = BATdescriptor(*yBAT_id)) == NULL) {
+ BBPreleaseref(xBAT->batCacheid);
+ throw(MAL, "batgeom.Distance", RUNTIME_OBJECT_MISSING);
+ }
+
+ //check if the BATs have dense heads and are aligned
+ if (!BAThdense(xBAT) || !BAThdense(yBAT)) {
+ ret = createException(MAL, "batgeom.Distance", "BATs must have
dense heads");
+ goto clean;
+ }
+ if(xBAT->hseqbase != yBAT->hseqbase || BATcount(xBAT) !=
BATcount(yBAT)) {
+ ret=createException(MAL, "batgeom.Distance", "BATs must be
aligned");
+ goto clean;
+ }
+
+ //get the GEOS representation of the geometry
+ if(!(geosGeometry = wkb2geos(*geomWKB)))
+ return createException(MAL, "batgeom.Contains", "wkb2geos
failed");
+
+ //chech the type of the geometry and choose the appropriate distance
function
+ switch(GEOSGeomTypeId(geosGeometry)+1) {
+ case wkbPoint:
+ if((outBAT = point2point_distance(geosGeometry, xBAT, yBAT)) ==
NULL) {
+ ret = createException(MAL, "batgeom.Distance", "Problem
evalauting the contains");
+ goto clean;
+ }
+ break;
+ default:
+ return createException(MAL, "batgeom.Distance", "This Geometry
type is not supported");
+ }
+
+ BBPkeepref(*outBAT_id = outBAT->batCacheid);
+ goto clean;
+
+clean:
+ if(xBAT)
+ BBPreleaseref(xBAT->batCacheid);
+ if(yBAT)
+ BBPreleaseref(yBAT->batCacheid);
+ return ret;
+}
+
/*
static BAT* BATMBRfilter(double xmin, double ymin, double xmax, double ymax,
wkb** geomWKB, int srid) {
BAT *outBAT = NULL, *leftBottomBAT = NULL, *leftTopBAT = NULL,
*rightBottomBAT = NULL, *rightTopBAT = NULL;
@@ -808,7 +905,7 @@ str wkbFilterWithImprints_geom_bat(bat*
return createException(MAL,"batgeom.Filter","Problem filtering
yBAT");
}
-fprintf(stderr, "Original MBR contains %u points\n", (unsigned
int)BATcount(candidateOIDsBAT));
+//fprintf(stderr, "Original MBR contains %u points\n", (unsigned
int)BATcount(candidateOIDsBAT));
//BATMBRfilter(xmin, ymin, xmax, ymax, geomWKB, (*geomWKB)->srid);
BBPreleaseref(xBAT->batCacheid);
BBPreleaseref(yBAT->batCacheid);
diff --git a/geom/sql/40_geom.sql b/geom/sql/40_geom.sql
--- a/geom/sql/40_geom.sql
+++ b/geom/sql/40_geom.sql
@@ -4058,7 +4058,8 @@ CREATE FUNCTION ST_Overlaps(geom1 Geomet
CREATE FUNCTION ST_Relate(geom1 Geometry, geom2 Geometry,
intersection_matrix_pattern string) RETURNS boolean EXTERNAL NAME geom."Relate";
--Distance between Geometries
CREATE FUNCTION ST_Distance(geom1 Geometry, geom2 Geometry) RETURNS double
EXTERNAL NAME geom."Distance";
-CREATE FUNCTION ST_Distance(geom1 Geometry, xCoordinate double, yCoordinate
double, srid integer) RETURNS double EXTERNAL NAME geom."Distance";
+CREATE FUNCTION ST_Distance1(geom1 Geometry, xCoordinate double, yCoordinate
double, srid integer) RETURNS double EXTERNAL NAME geom."Distance1";
+CREATE FUNCTION ST_Distance2(geom1 Geometry, xCoordinate double, yCoordinate
double, srid integer) RETURNS double EXTERNAL NAME geom."Distance2";
--Functions that implement spatial operators
CREATE FUNCTION ST_Intersection(geom1 Geometry, geom2 Geometry) RETURNS
Geometry EXTERNAL NAME geom."Intersection";
CREATE FUNCTION ST_Difference(geom1 Geometry, geom2 Geometry) RETURNS Geometry
EXTERNAL NAME geom."Differnce";
diff --git a/monetdb5/optimizer/opt_geospatial.c
b/monetdb5/optimizer/opt_geospatial.c
--- a/monetdb5/optimizer/opt_geospatial.c
+++ b/monetdb5/optimizer/opt_geospatial.c
@@ -118,8 +118,12 @@ int OPTgeospatialImplementation(Client c
actions +=2;
}
- } else if(strcasecmp(getFunctionId(oldInstrPtr[i]),
"distance") == 0 && strcasecmp(getFunctionId(oldInstrPtr[i+1]),
"thetasubselect") == 0) {
- //I should check the theta comparison. In case
it is > OR >= then there should be no filtering
+ } else if((strcasecmp(getFunctionId(oldInstrPtr[i]),
"distance1") == 0 || strcasecmp(getFunctionId(oldInstrPtr[i]), "distance2") ==
0)
+ &&
strcasecmp(getFunctionId(oldInstrPtr[i+1]), "thetasubselect") == 0) {
+ //the filter does not make sense if comparison
is > OR >=
+ // if(strcmp(getArg(oldInstrPtr[i+1],4), ">")==0
|| strcmp(getArg(oldInstrPtr[i+1],4),">=")==0) {
+ // pushInstruction(mb, oldInstrPtr[i]);
+ // } else
if(oldInstrPtr[i]->argc == 5) {
InstrPtr bufferInstrPtr;
int bufferReturnId;
_______________________________________________
checkin-list mailing list
[email protected]
https://www.monetdb.org/mailman/listinfo/checkin-list
