Changeset: eea7c545ffe7 for MonetDB
URL: https://dev.monetdb.org/hg/MonetDB/rev/eea7c545ffe7
Modified Files:
geom/monetdb5/geom.c
geom/monetdb5/geom.h
Branch: geo-update
Log Message:
Refactored code, improved Point-In-Polygon calculation
diffs (truncated from 1066 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
@@ -15,28 +15,33 @@
#include "gdk_logger.h"
#include "mal_exception.h"
+static wkb *geos2wkb(const GEOSGeometry *geosGeometry);
+
+/**
+*
+* Geographic update code start
+*
+**/
+
/**
- * Convertions
- *
- **/
+* Convertions
+*
+**/
const double earth_radius = 6371.009;
const double earth_radius_meters = 6371009;
-/**
- * Converts a longitude value in degrees to radians
- * The normalization part was taken from PostGIS
- */
+/* Converts a longitude value in degrees to radians */
static double deg2RadLongitude(double lon_degrees)
{
//Convert
double lon = M_PI * lon_degrees / 180.0;
//Normalize
+ //TODO PostGIS code, refactor
if (lon == -1.0 * M_PI)
return M_PI;
if (lon == -2.0 * M_PI)
return 0.0;
-
if (lon > 2.0 * M_PI)
lon = remainder(lon, 2.0 * M_PI);
@@ -55,15 +60,13 @@ static double deg2RadLongitude(double lo
return lon;
}
-/**
- * Converts a latitude value in degrees to radians
- * The normalization part was taken from PostGIS
- */
+/* Converts a latitude value in degrees to radians */
static double deg2RadLatitude(double lat_degrees)
{
//Convert
double lat = M_PI * lat_degrees / 180.0;
//Normalize
+ //TODO PostGIS code, refactor
if (lat > 2.0 * M_PI)
lat = remainder(lat, 2.0 * M_PI);
@@ -93,17 +96,15 @@ static GeoPoint deg2RadPoint(GeoPoint ge
return geo;
}
-#if 0
-
/**
* Converts a longitude value in radians to degrees
- * The normalization part was taken from PostGIS
*/
static double rad2DegLongitude(double lon_radians)
{
//Convert
double lon = lon_radians * 180.0 / M_PI;
//Normalize
+ //TODO PostGIS code, refactor
if (lon > 360.0)
lon = remainder(lon, 360.0);
@@ -127,13 +128,13 @@ static double rad2DegLongitude(double lo
/**
* Converts a latitude value in radians to degrees
- * The normalization part was taken from PostGIS
*/
static double rad2DegLatitude(double lat_radians)
{
//Convert
double lat = lat_radians * 180.0 / M_PI;
//Normalize
+ //TODO PostGIS code, refactor
if (lat > 360.0)
lat = remainder(lat, 360.0);
@@ -154,14 +155,16 @@ static double rad2DegLatitude(double lat
return lat;
}
+
/* Converts the GeoPoint from degrees to radians latitude and longitude*/
-static void rad2DegPoint(GeoPoint *geo)
-{
- geo->lon = rad2DegLongitude(geo->lon);
- geo->lat = rad2DegLatitude(geo->lat);
-}
-#endif
-
+static GeoPoint rad2DegPoint(GeoPoint geo)
+{
+ geo.lon = rad2DegLongitude(geo.lon);
+ geo.lat = rad2DegLatitude(geo.lat);
+ return geo;
+}
+
+/* Converts the a GEOSGeom Point into a GeoPoint */
static GeoPoint geoPointFromGeom(GEOSGeom geom)
{
GeoPoint geo;
@@ -170,16 +173,7 @@ static GeoPoint geoPointFromGeom(GEOSGeo
return geo;
}
-/*static GeoLine geoLineFromGeom(GEOSGeom geom)
-{
- GeoLine geo;
- GEOSGeomGetX(GEOSGeomGetStartPoint(geom), &(geo.start.lon));
- GEOSGeomGetY(GEOSGeomGetStartPoint(geom), &(geo.start.lat));
- GEOSGeomGetX(GEOSGeomGetEndPoint(geom), &(geo.end.lon));
- GEOSGeomGetY(GEOSGeomGetEndPoint(geom), &(geo.end.lat));
- return geo;
-}*/
-
+/* Converts the a GEOSGeom Line into GeoLines (one or more line segments) */
static GeoLines geoLinesFromGeom(GEOSGeom geom)
{
const GEOSCoordSequence *gcs = GEOSGeom_getCoordSeq(geom);
@@ -193,6 +187,7 @@ 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);
}
+ geo.bbox = NULL;
return geo;
}
@@ -230,7 +225,7 @@ static str wkbGetComplatibleGeometries(w
{
err = createException(MAL, "geom.wkbGetComplatibleGeometries",
SQLSTATE(38000) "Geos operation wkb2geos failed");
}
- else if (GEOSGetSRID((*ga)) != GEOSGetSRID((*gb)))
+ else if (GEOSGetSRID((*ga)) != GEOSGetSRID(*gb))
{
GEOSGeom_destroy(*ga);
GEOSGeom_destroy(*gb);
@@ -240,7 +235,8 @@ static str wkbGetComplatibleGeometries(w
}
/**
-* Convert spherical coordinates to cartesian coordinates on unit sphere
+* Convert spherical coordinates to cartesian coordinates on unit sphere.
+* The inputs have to be in radians.
*/
static CartPoint geo2cart(GeoPoint geo)
{
@@ -252,20 +248,15 @@ static CartPoint geo2cart(GeoPoint geo)
}
/**
-* Convert spherical coordinates to cartesian coordinates
+* Convert spherical coordinates to cartesian coordinates on unit sphere.
+* The inputs have to be in degrees.
*/
-static CartPoint geo2cart_r(GeoPoint geo)
-{
- CartPoint cart;
- cart.x = earth_radius * cos(geo.lat) * cos(geo.lon);
- cart.y = earth_radius * cos(geo.lat) * sin(geo.lon);
- cart.z = earth_radius * sin(geo.lat);
- return cart;
-}
-
-/**
-* Convert cartesian coordinates to spherical coordinates on unit sphere
-*/
+static CartPoint geo2cartFromDegrees(GeoPoint geo)
+{
+ return geo2cart(deg2RadPoint(geo));
+}
+
+/* Convert cartesian coordinates to spherical coordinates on unit sphere */
static GeoPoint cart2geo(CartPoint cart)
{
GeoPoint geo;
@@ -274,25 +265,163 @@ static GeoPoint cart2geo(CartPoint cart)
return geo;
}
+/* Converts two lat/lon points into cartesian coordinates and creates a Line
geometry */
+static GEOSGeom cartesianLineFromGeoPoints(GeoPoint p1, GeoPoint p2)
+{
+ CartPoint p1_cart, p2_cart;
+ p1_cart = geo2cartFromDegrees(p1);
+ p2_cart = geo2cartFromDegrees(p2);
+ GEOSCoordSequence *lineSeq = GEOSCoordSeq_create(2, 3);
+ GEOSCoordSeq_setXYZ(lineSeq, 0, p1_cart.x, p1_cart.y, p1_cart.z);
+ GEOSCoordSeq_setXYZ(lineSeq, 1, p2_cart.x, p2_cart.y, p2_cart.z);
+ return GEOSGeom_createLineString(lineSeq);
+}
+
/**
-* Convert cartesian coordinates to spherical coordinates
-*/
-/*static GeoPoint cart2geo_r(CartPoint cart)
-{
- GeoPoint geo;
- geo.lon = atan2(cart.y, cart.x);
- geo.lat = asin(cart.z / earth_radius);
- return geo;
-}*/
+* Bounding Box functions
+*
+**/
+/* Adds a Cartesian Point to the BoundingBox */
+static void boundingBoxAddPoint(BoundingBox *bb, CartPoint p)
+{
+ if (bb->xmin > p.x)
+ bb->xmin = p.x;
+ if (bb->xmax < p.x)
+ bb->xmax = p.x;
+ if (bb->ymin > p.y)
+ bb->ymin = p.y;
+ if (bb->ymax < p.y)
+ bb->ymax = p.y;
+ if (bb->zmin > p.z)
+ bb->zmin = p.z;
+ if (bb->zmax < p.z)
+ bb->zmax = p.z;
+}
+
+/* Builds the BoundingBox for a GeoLines geometry */
+static BoundingBox* boundingBoxLines(GeoLines lines)
+{
+ CartPoint c;
+ BoundingBox *bb = GDKzalloc(sizeof(BoundingBox));
+
+ //If there are no segments, return NULL
+ if (lines.segmentCount == 0)
+ {
+ return NULL;
+ }
+
+ c = geo2cartFromDegrees(lines.segments[0].start);
+
+ //Initialize the bounding box with the first point
+ bb->xmin = bb->xmax = c.x;
+ bb->ymin = bb->ymax = c.y;
+ bb->zmin = bb->zmax = c.z;
+
+ for (int i = 0; i < lines.segmentCount; i++)
+ {
+ c = geo2cartFromDegrees(lines.segments[i].end);
+ boundingBoxAddPoint(bb, c);
+ }
+ return bb;
+}
+
+static int boundingBoxContainsPoint(BoundingBox bb, CartPoint pt)
+{
+ 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));
+ copy->xmin = bb.xmin;
+ copy->xmax = bb.xmax;
+ copy->ymin = bb.ymin;
+ copy->ymax = bb.ymax;
+ copy->zmin = bb.zmin;
+ copy->zmax = bb.zmax;
+ return *copy;
+}
+
+/* Returns a point outside of the polygon's bounding box, for Point-In-Polygon
calculation */
+static GeoPoint pointOutsidePolygon(GeoLines polygonRing)
+{
+ //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: From POSTGIS -> CHANGE
+ double grow = M_PI / 180.0 / 60.0;
+ CartPoint corners[8];
+ while (grow < M_PI)
+ {
+ if (bb.xmin > -1)
_______________________________________________
checkin-list mailing list
[email protected]
https://www.monetdb.org/mailman/listinfo/checkin-list
