Changeset: 79f196ea8e2c for MonetDB
URL: https://dev.monetdb.org/hg/MonetDB/rev/79f196ea8e2c
Modified Files:
geom/monetdb5/geom.c
geom/monetdb5/geom.h
Branch: geo-update
Log Message:
Added GeoPolygon data type to also represent the interior rings of a Polygon.
Added GeoPolygon to all Polygon functions in geom.c, corrected the
pointInPolygon error of false positives if the point was within a polygon hole.
diffs (truncated from 516 to 300 lines):
diff --git a/geom/monetdb5/geom.c b/geom/monetdb5/geom.c
--- a/geom/monetdb5/geom.c
+++ b/geom/monetdb5/geom.c
@@ -174,7 +174,7 @@ static GeoPoint geoPointFromGeom(GEOSGeo
}
/* Converts the a GEOSGeom Line into GeoLines (one or more line segments)
- Argument must be a Line, use geoLinesFromGeomPolygon for Polygons. */
+ Argument must be a Line geometry. */
static GeoLines geoLinesFromGeom(GEOSGeom geom)
{
const GEOSCoordSequence *gcs = GEOSGeom_getCoordSeq(geom);
@@ -188,13 +188,33 @@ static GeoLines geoLinesFromGeom(GEOSGeo
GEOSCoordSeq_getXY(gcs, i, &geo.segments[i].start.lon,
&geo.segments[i].start.lat);
GEOSCoordSeq_getXY(gcs, i + 1, &geo.segments[i].end.lon,
&geo.segments[i].end.lat);
}
+ //TODO Calculate Boundind Box on initializion?
geo.bbox = NULL;
return geo;
}
-static GeoLines geoLinesFromGeomPolygon(GEOSGeom geom)
-{
- return geoLinesFromGeom((GEOSGeom)GEOSGetExteriorRing(geom));
+/* Converts the a GEOSGeom Line into GeoPolygon (with exterior ring and
zero-to-multiple interior rings)
+ Argument must be a Polygon geometry. */
+static GeoPolygon geoPolygonFromGeom(GEOSGeom geom)
+{
+ GeoPolygon geo;
+ //TODO Calculate Boundind Box on initializion?
+ geo.bbox = NULL;
+ //Get exterior ring GeoLines
+ geo.exteriorRing =
geoLinesFromGeom((GEOSGeom)GEOSGetExteriorRing(geom));
+ geo.interiorRingsCount = GEOSGetNumInteriorRings(geom);
+ //If there are interior rings, allocate space to their GeoLines
representation
+ if (geo.interiorRingsCount > 0)
+ {
+ //TODO Malloc fail exception?
+ geo.interiorRings = GDKmalloc(sizeof(GeoLines) *
geo.interiorRingsCount);
+ }
+ //Get interior rings GeoLines
+ for (int i = 0; i < geo.interiorRingsCount; i++)
+ {
+ geo.interiorRings[i] =
geoLinesFromGeom((GEOSGeom)GEOSGetInteriorRingN(geom, i));
+ }
+ return geo;
}
static GeoPoint geoPointFromLatLon(double lon, double lat)
@@ -305,7 +325,7 @@ static void boundingBoxAddPoint(Bounding
}
/* Builds the BoundingBox for a GeoLines geometry */
-static BoundingBox* boundingBoxLines(GeoLines lines)
+static BoundingBox *boundingBoxLines(GeoLines lines)
{
CartPoint c;
BoundingBox *bb = GDKzalloc(sizeof(BoundingBox));
@@ -336,8 +356,9 @@ static int boundingBoxContainsPoint(Boun
return bb.xmin <= pt.x && bb.xmax >= pt.x && bb.ymin <= pt.y && bb.ymax
>= pt.y && bb.zmin <= pt.z && bb.zmax >= pt.z;
}
-static BoundingBox boundingBoxCopy (BoundingBox bb) {
- BoundingBox* copy = GDKmalloc(sizeof(BoundingBox));
+static BoundingBox boundingBoxCopy(BoundingBox bb)
+{
+ BoundingBox *copy = GDKmalloc(sizeof(BoundingBox));
copy->xmin = bb.xmin;
copy->xmax = bb.xmax;
copy->ymin = bb.ymin;
@@ -348,14 +369,16 @@ static BoundingBox boundingBoxCopy (Boun
}
/* Returns a point outside of the polygon's bounding box, for Point-In-Polygon
calculation */
-static GeoPoint pointOutsidePolygon(GeoLines polygonRing)
+static GeoPoint pointOutsidePolygon(GeoPolygon polygon)
{
//If the geometry doesn't have its BoundingBox calculated, calculate it
- if (polygonRing.bbox == NULL) {
- polygonRing.bbox = boundingBoxLines(polygonRing);
- }
- BoundingBox bb = *polygonRing.bbox;
- BoundingBox bb2 = boundingBoxCopy(*polygonRing.bbox);
+ //TODO Should we consider the interior rings for bounding box
calculation?
+ if (polygon.bbox == NULL)
+ {
+ polygon.bbox = boundingBoxLines(polygon.exteriorRing);
+ }
+ BoundingBox bb = *polygon.bbox;
+ BoundingBox bb2 = boundingBoxCopy(*polygon.bbox);
//TODO: From POSTGIS -> CHANGE
double grow = M_PI / 180.0 / 60.0;
@@ -417,6 +440,7 @@ static GeoPoint pointOutsidePolygon(GeoL
}
grow *= 2.0;
}
+ //TODO: Should this be the return value in case no point is found?
return geoPointFromLatLon(0, 0);
}
@@ -545,53 +569,39 @@ static double geoDistanceLineLine(GeoLin
}
//TODO Implement intersection ourselves so we don't use GEOS?
-//TODO This still produces some false positives -> Is it only for Polygons
where the Point is in a hole?
/* Checks if a Point is within a Polygon */
-static bool pointWithinPolygon(GeoLines polygonRing,GeoPoint point)
+static bool pointWithinPolygon(GeoPolygon polygon, GeoPoint point)
{
int intersectionNum = 0;
GEOSGeometry *segmentPolygon, *intersectionPoints;
+ GeoLines polygonRing;
//Get an point that's outside the polygon
- GeoPoint outsidePoint = pointOutsidePolygon(polygonRing);
+ GeoPoint outsidePoint = pointOutsidePolygon(polygon);
//No outside point was found, return false
- if (outsidePoint.lat == 0 && outsidePoint.lon == 0) {
+ if (outsidePoint.lat == 0 && outsidePoint.lon == 0)
+ {
return false;
}
-
+
/*printf("Outside point: (%f %f)\n",outsidePoint.lon, outsidePoint.lat);
fflush(stdout);*/
//Construct a line between the outside point and the input point
GEOSGeometry *outInLine = cartesianLineFromGeoPoints(point,
outsidePoint);
- //Count the number of intersections between the polygon and the
constructed line
+ //Count the number of intersections between the polygon exterior ring
and the constructed line
+ polygonRing = polygon.exteriorRing;
for (int i = 0; i < polygonRing.segmentCount; i++)
{
-
segmentPolygon =
cartesianLineFromGeoPoints(polygonRing.segments[i].start,
polygonRing.segments[i].end);
intersectionPoints = GEOSIntersection(segmentPolygon,
outInLine);
- //printf("Segment (%d): (%f %f)->(%f
%f)\n",i,polygonRing.segments[i].start.lon, polygonRing.segments[i].start.lat,
polygonRing.segments[i].end.lon, polygonRing.segments[i].end.lat);
//If there is an intersection, a point will be returned (line
when there is none)
if (GEOSGeomTypeId(intersectionPoints) == GEOS_POINT)
{
intersectionNum++;
- /*CartPoint pDegrees;
- GeoPoint pRadians;
- double x, y, z;
- GEOSGeomGetX(intersectionPoints, &x);
- GEOSGeomGetY(intersectionPoints, &y);
- GEOSGeomGetZ(intersectionPoints, &z);
- pDegrees.x = x;
- pDegrees.y = y;
- pDegrees.z = z;
- pRadians = rad2DegPoint(cart2geo(pDegrees));
- printf("Intersection Num %d on Segment Num %d\n",
intersectionNum,i);
- printf("Intersection Point (Degrees) (%f %f)\n",
pRadians.lon, pRadians.lat);
- printf("Line (%f %f)->(%f %f)\nLine (%f %f)->(%f
%f)\n", outsidePoint.lon, outsidePoint.lat, point.lon, point.lat,
polygonRing.segments[i].start.lon, polygonRing.segments[i].start.lat,
polygonRing.segments[i].end.lon, polygonRing.segments[i].end.lat);
- fflush(stdout);*/
}
if (intersectionPoints != NULL)
@@ -599,6 +609,29 @@ static bool pointWithinPolygon(GeoLines
if (segmentPolygon != NULL)
GEOSGeom_destroy(segmentPolygon);
}
+
+ //Count the number of intersections between the polygon interior rings
and the constructed line
+ for (int j = 0; j < polygon.interiorRingsCount; j++)
+ {
+ polygonRing = polygon.interiorRings[j];
+ for (int i = 0; i < polygonRing.segmentCount; i++)
+ {
+ segmentPolygon =
cartesianLineFromGeoPoints(polygonRing.segments[i].start,
polygonRing.segments[i].end);
+ intersectionPoints = GEOSIntersection(segmentPolygon,
outInLine);
+
+ //If there is an intersection, a point will be returned
(line when there is none)
+ if (GEOSGeomTypeId(intersectionPoints) == GEOS_POINT)
+ {
+ intersectionNum++;
+ }
+
+ if (intersectionPoints != NULL)
+ GEOSGeom_destroy(intersectionPoints);
+ if (segmentPolygon != NULL)
+ GEOSGeom_destroy(segmentPolygon);
+ }
+ }
+
if (outInLine != NULL)
GEOSGeom_destroy(outInLine);
@@ -607,27 +640,35 @@ static bool pointWithinPolygon(GeoLines
}
/* Distance between Point and Polygon.*/
-static double geoDistancePointPolygon(GeoPoint point, GeoLines polygonRing)
+static double geoDistancePointPolygon(GeoPoint point, GeoPolygon polygon)
{
//Check if point is in polygon
- if (pointWithinPolygon(polygonRing,point))
- {
+ if (pointWithinPolygon(polygon, point))
return 0;
- }
-
- //Compare Point to the various polygon segments
- return geoDistancePointLine(point, polygonRing);
+
+ //Calculate distance from Point to the exterior and interior rings of
the polygon
+ double distance, min_distance = INT_MAX;
+ //First, calculate distance to the exterior ring
+ min_distance = geoDistancePointLine(point, polygon.exteriorRing);
+ //Then, calculate distance to the interior rings
+ for (int i = 0; i < polygon.interiorRingsCount; i++)
+ {
+ distance = geoDistancePointLine(point,
polygon.interiorRings[i]);
+ if (distance < min_distance)
+ min_distance = distance;
+ }
+ return min_distance;
}
/* Distance between Line and Polygon. */
-static double geoDistanceLinePolygon(GeoLines line, GeoLines polygon)
+static double geoDistanceLinePolygon(GeoLines line, GeoPolygon polygon)
{
double distance, min_distance = INT_MAX;
//Calculate distance to all start vertices of the line
for (int i = 0; i < line.segmentCount; i++)
{
distance = geoDistancePointPolygon(line.segments[i].start,
polygon);
-
+
//Short-cut in case the point is within the polygon
if (distance == 0)
return 0;
@@ -637,19 +678,18 @@ static double geoDistanceLinePolygon(Geo
}
//Calculate distance to the last vertice (not covered by the previous
loop)
distance = geoDistancePointPolygon(line.segments[line.segmentCount -
1].end, polygon);
- if (distance < min_distance)
- min_distance = distance;
- return min_distance;
+ return distance < min_distance ? distance : min_distance;
}
/* Distance between two Polygons. */
-static double geoDistancePolygonPolygon(GeoLines polygon1, GeoLines polygon2)
+//TODO Does this calculate the minimum distance between the interior rings as
well? We are comparing the exterior ring of each with all segments of the other.
+static double geoDistancePolygonPolygon(GeoPolygon polygon1, GeoPolygon
polygon2)
{
double distance1, distance2;
- //Calculate the distance between all vertices of polygon1 and segments
of polygon2
- distance1 = geoDistanceLinePolygon(polygon1, polygon2);
- //Calculate the distance between all vertices of polygon2 and segments
of polygon1
- distance2 = geoDistanceLinePolygon(polygon2, polygon1);
+ //Calculate the distance between the exterior ring of polygon1 and all
segments of polygon2 (including the interior rings)
+ distance1 = geoDistanceLinePolygon(polygon1.exteriorRing, polygon2);
+ //Other way around
+ distance2 = geoDistanceLinePolygon(polygon2.exteriorRing, polygon1);
//And return the minimum
return distance1 < distance2 ? distance1 : distance2;
}
@@ -683,27 +723,27 @@ static double geoDistanceSingle(GEOSGeom
else if (dimA == 0 && dimB == 2)
{
/* Point and Polygon */
- return geoDistancePointPolygon(geoPointFromGeom(a),
geoLinesFromGeomPolygon(b));
+ return geoDistancePointPolygon(geoPointFromGeom(a),
geoPolygonFromGeom(b));
}
else if (dimA == 2 && dimB == 0)
{
/* Polygon and Point */
- return geoDistancePointPolygon(geoPointFromGeom(b),
geoLinesFromGeomPolygon(a));
+ return geoDistancePointPolygon(geoPointFromGeom(b),
geoPolygonFromGeom(a));
}
else if (dimA == 1 && dimB == 2)
{
/* Line/LinearRing and Polygon */
- return geoDistanceLinePolygon(geoLinesFromGeom(a),
geoLinesFromGeomPolygon(b));
+ return geoDistanceLinePolygon(geoLinesFromGeom(a),
geoPolygonFromGeom(b));
}
else if (dimA == 2 && dimB == 1)
{
/* Polygon and Line/LinearRing */
- return geoDistanceLinePolygon(geoLinesFromGeom(b),
geoLinesFromGeomPolygon(a));
+ return geoDistanceLinePolygon(geoLinesFromGeom(b),
geoPolygonFromGeom(a));
}
else if (dimA == 2 && dimB == 2)
{
/* Polygon and Polygon */
- return geoDistancePolygonPolygon(geoLinesFromGeomPolygon(a),
geoLinesFromGeomPolygon(b));
+ return geoDistancePolygonPolygon(geoPolygonFromGeom(a),
geoPolygonFromGeom(b));
}
return INT_MAX;
}
@@ -793,7 +833,8 @@ str wkbIntersectsGeographic(bit *out, wk
}
_______________________________________________
checkin-list mailing list
[email protected]
https://www.monetdb.org/mailman/listinfo/checkin-list
